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PREFACE 

Calculation procedures for nonreacting compressible two- and three- 
dimensional turbulent boundary layers are reviewed. A summary of integral, 
transformation, and correlation methods, as well as finite-difference solutions 
of the complete boundary-layer equations is included. Alternative numerical 
solution procedures are examined, and both mean field and mean turbulence field 
closure models are considered. A discussion of physics and related calculation 
problems peculiar to compressible turbulent boundary layers is included. A 
listing of available solution procedures (finite-difference, finite-element, and 
weighted-residual methods) is provided. Detailed consideration is given to 
influence of compressibility, low Reynolds number, wall blowing, and pressure 
gradient upon mean field closure constants. 

The information contained in this publication was presented at the 1976 
Von Karman Institute for Fluid Dynamics Lecture Series "Compressible Turbulent 
Boundary Layers," March 1-5, 1976, Rhode-St.-Genese , Belgium. 
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INTRODUCTION 


The design of missiles, manned and unmanned entry vehicles, and transporta- 
tion systems capable of speeds in the transonic and supersonic regimes requires 
detailed information on such boundary- layer quantities as skin friction, aero- 
dynamic heating, and viscous displacement thickness and mass flow. Typical com- 
ponents requiring such design information include nacelles, control surfaces, 
turbomachinery blading, nozzles, airfoils and fuselage, inlets, and combustors. 
Detailed boundary-layer predictions are also needed for the design of facilities 
used in experimental investigations. 

The various stages of design demand boundary-layer information in increas- 
ing detail. For this reason and also because flow complexity varies considera- 
bly (depending upon the design and component) , a heirarchy of calculation pro- 
cedures has been developed over the years, ranging from simple, "back of the 
envelope" methods to complex, numerical approaches which require use of the 
largest digital computers. In the present paper this entire gamut of methods is 
reviewed, with emphasis on the more complete procedures which solve numerically 
the partial-differential equations governing boundary-layer motion and on the 
influence of conditions usually encountered in practice (such as pressure gradi- 
ent, mass injection, and low Reynolds number) upon the necessary closure "con- 
stants" used in representing the turbulent shear. Also included is a detailed 
discussion of 30 me physics and resultant calculation problems which are peculiar 
to compressible turbulent boundary layers. The basic purpose of this volume is 
to provide a ready reference and introduction to the various procedures cur- 
rently available for calculation of compressible turbulent boundary layers. 
Therefore, also included is a listing of available methods of the more complete 
type and some discussion of the various alternate numerical procedures which can 
be used for solving the nonlinear partial-differential equations governing fluid 
motion in compressible turbulent boundary layers. This review does not include 
detailed consideration of time-dependent boundary layers, relaminarization, and 
heterogeneous or chemically reacting flows. 


SYMBOLS 

A Van Driest damping factor (eq. (41)) 

A + = Au x /v 

a speed of sound 

a-) = u'v'/2e 

Cf skin friction coefficient, t w ^/1 p e u e 2 

Cf >0 skin friction coefficient without wall blowing 


c 


wing chord 




%u Nusselt number, (Ngt) (Npr )(^e ,x) 

Np r molecular Prandtl number 

Np r> T total turbulent Prandtl number, e/<x 
Np r ,t static turbulent Prandtl number, e/<t 

Ngt Stanton number, q w /p e u e c p (T aw - T w ) 

p pressure 

p + pressure gradient parameter, (vu e /u x 3) (du e /dx) 

q wall heat transfer rate 


2 


R 

r c 

r N 

T 


u,v,w 


u i 

u + 

u« 

U T 

V 

x,y,z 

x i 

y + 

a 

r 

r x 

r y 

Y 
6 


6 * 


6+ 

e 

C 

Tl 


Reynolds number, pu/y; also universal gas constant 
transverse radius of curvature 
nose radius 
temperature 

velocity components in x-, y-, and z-directions 
general velocity notation (i = 1, 2, 3) 

= u/u T 

Van Driest* s generalized velocity 
friction velocity, (^w/Pw) 1 ^ 2 
total velocity vector 
Cartesian coordinates 
curvilinear coordinates (i = 1, 2, 3) 

= yu T /V 

Clauser constant; also angle of attack 

* 

pressure gradient parameter, (6 ^/t w ) ( dp/dx) 

intermittency function 

streamwise intermittency function 

normal intermittency function 

ratio of specific heats 

boundary-layer thickness 

1 - P u \ dy 
Pe u e/ 

= <Su T /v w 

dynamic eddy viscosity 

transformed normal coordinate, Crocco (eq. (91)) 

transformed normal coordinate, Levy-Lees (eq. (32)); also transformed 
transverse coordinate for Crocco variables (see fig. 53) 


displacement thickness 


■ 

J n 


3 


0 


0 


k t 


K t 


= T/T e 


momentum thickness, 


_ v4i’ 
3H/8y 



dy 


_ v'h' 
8h/3y 


y molecular viscosity 

v kinematic viscosity 

£ transformed streamwise coordinate (eq. (31)); also x-| -coordinate for 

Crocco variables 

p density 

a kinematic eddy viscosity (eq. (40)) 

x shear stress 

$ shear function (eq. (92)) 

4> spreading angle (fig. 5) 

q> circumferential angle 

stream function 
co vorticity 

Subscripts: 

aw adiabatic wall 

bl boundary layer 

cr critical value 

e local edge of boundary layer 

eff effective value 

i incompressible 

i,j,k indices 

L reference length 
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max maximum value 

n index integer 

o stagnation value 

p region of peak v' 2 

t total 

tr transition value 

vo based on virtual origin location 

x based on distance coordinate x 

x^ coordinate direction (i = 1, 2, 3) 

w wall value 

6 based on momentum thickness 

00 free stream 

A prime denotes a fluctuating quantity. 

A bar indicates a mean quantity. 


MAJOR GENERAL REFERENCES 

The major reference books pertinent to the calculation of compressible tur- 
bulent boundary layers were mainly written in two distinct time frames: late 

fifties to early sixties (refs. 1 to 5) and late sixties to early seventies 
(refs. 6 to 12). Of these works, the ones most useful in the preparation of the 
present volume on calculation methods included references 2, 5, 8, 9, and 12. 
Several review and background articles are also available (e.g., refs. 13 to 25) 
All of these were quite valuable, especially the papers of Reynolds (refs. 19 
and 20) and Bradshaw (ref. 22). Another category of general references is con- 
ference proceedings (e.g., refs. 26 to 30). These are excellent sources, partic 
ularly for comparisons between data and theory (especially refs. 28 and 30). 
Reference 26 contains many of the fundamental concepts, such as the Morkovin 
hypothesis, which are the foundation for the current generation of compressible 
calculation procedures. A final category of general references includes reviews 
of available data (e.g., refs. 11 and 31 to 37), which are especially important 
for evaluation of test cases suitable to "calibrate" the various turbulence clos 
ure constants. 
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EQUATIONS AND CLOSURE METHODS 


Governing Differential Equations 

The basic differential equations governing compressible turbulent boundary- 
layer flow are (1) a statement of the conservation of mass, (2) momentum equa- 
tions obtained from the Navier-Stokes equations, and (3) an expression for the 
conservation of energy. Also needed for solution is an equation of state (ideal 
gas assumed herein) and equations for molecular transport properties. 

Given this set of equations it is conceptually possible to integrate for- 
ward in time on a computer and, given sufficient grid resolution, obtain the tur- 
bulent motions "exactly" with very little empiricism. However, as noted in ref- 
erence 12 and elsewhere, the computer capacity currently available falls several 
orders of magnitude short of the capability needed for such a solution. There- 
fore, one must resort to the usual Reynolds averaging, where the flow is repre- 
sented by a time mean and an instantaneous fluctuation; for example, 

u=u+u' v=v+v* 

where the following averaging rules apply: 


u v* = 0 u'v' ^ 0 

To simplify the problem for ease of solution, with very little loss of accuracy 
(ref. 38), order of magnitude arguments are made (usually referred to as the 
boundary-layer assumptions); that is, 


x >> y u » v 


L_ « 9_ 

3x 3y 


The resultant set of governing equations in surface curvilinear coordinates for 
the steady, two-dimensional boundary-layer flow (6/r c « 1) of a compressible 
ideal gas (in the turbulent state) are 


Conservation of mass: 


3(pu) + 3(pv) = 0 
9x 3y 

Longitudinal momentum: 

p 3u _ (pv) 'u' 

3y 


pu + pv iLy = - ^2 + 

3x 3y 3x 3y 


( 1 ) 


( 2 ) 


Normal momentum: 



Conservation of energy: 


pu M + pv M = L 
9x 3y 3y 



3H 

_9y 


+ 


(N Pr - 


1)9u£_ 
2 3y 


(pv)'h' 



(4) 


Detailed discussion and derivation of these equations are available in refer- 
ences 2, 12, 39, 40, and 41 . The equations are also given in references 42 and 
43, in terms of mass-weighted dependent variables. 


There are two points to be made concerning these basic equations (eqs. (1) 
to (4)). First of all, in the averaging process the usual new unknowns have 
appeared (shown underlined) which account for the turbulent fluxes of momentum 
and energy. The specification of these quantities, in terms of known parame- 
ters, comprises the major difficulty in the calculation of compressible turbu- 
lent boundary layers and is usually termed the "closure problem." The other 
point concerns the term pv. Expanded out, this term is 

pv = p v + p’v' (5) 


Since p'v' appears in the equations with normal or y-deriv itive s, the 
term is not negligible and must be accounted for. Fortunately, p'v' always 
appears with p v, and therefore, a new definition of v can be used 

v = v+P5! (6) 

P 


(which is actually a mass-weighted variable), and the influence of p'v' can 
be included implicitly in the solution. This approach is satisfactory as long 
as the act ual v value is not required. If v values must be computed, some 
model of p'v* is obviously necessary. This inclusion of p'v* is an impor- 
tant issue. If the pv term were not handled in the manner shown in equ a- 
tions (2) and (4), then two Reynolds stre ss t er ms ap pear: u p'v' and p u'v'. 
Morkovin, in reference 44, shows that u p'v’ /p u'v' (or the ratio of the com- 
pressibility term to the usual low-speed Reynolds stress term) can b e of the 
order of 0.6 to 1.0 for a supersonic boundary layer, and therefore p'v' would 
have to be accurately modeled (an extremely difficult task). 


Closure Methods 

In the present paper the usual breakdown of closure procedures (e.g. , 
ref. 45) into (1) simple or zeroth-order methods, (2) first-order or mean field 
closure methods, and (3) second-order or mean turbulence field closure methods 
is followed . 

The zeroth-order case consists of two major subcategories: integral 

approaches and empirical laws, such as Cf correlations. In the most general 
of these approaches, the integral methods, equations (1) to (4) are formally 
integrated in the normal or y-direction. As a result of this procedure the 
unknown turbulent flux terms disappear, but their influence is still present in 
that profiles must be supplied (assumed, obtained from data, etc.), and these pro- 
files are influenced to a great extent by the turbulence induced fluxes. There- 
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fore , the simple or zeroth-order procedures are characterized by a requirement 
for substantial amounts of empiricism. 

In the mean field closure approaches the partial-differential equations are 
solved directly (eqs. (1) to (4)) and the turbulent flux terms are related to 
mean flow quantities. This approach is nearly correct (and indeed quite exact) 
for slowly varying flows and over a wide variety of boundary conditions. At the 
present time this closure approach is tending to supplant the integral approaches 
in industrial application, primarily as a result of a greater reliability (accu- 
racy of prediction) over a wide range of conditions and routine use of high- 
speed computers. This increased capability is purchased at the expense of con- 
siderably increased computer time (compared with the best of the integral 
methods) . 

In mean turbulent field closure, the differential equations (derived from 
the Navier-Stokes equations, as described in a later section) governing the tur- 
bulence flux terms are solved. These new equations involve additional unknowns, 
but the mean profiles generally are relatively insensitive to the precise 
details of the modeling in these second-order equations. This approach has actu- 
ally had limited application to the compressible boundary-layer case (primarily 
as a result of the success of the mean field methods) and is really needed only 
in cases where the flow undergoes a sudden change in boundary condition or expe- 
riences a large gradient (a situation generally termed nonequilibrium). 

The usual range of application of each of these three closure approaches is 
indicated schematically as follows: 


EQUILIBRIUM 


NEAR EQUILIBRIUM 


NONEQUILIBRIUM 


•ZERO OR EQUILIBRIUM 
GRADIENTS OF 
PRESSURE. WALL 
TEMPERATURE. AND 
WALL INJECTION 


•MODERATE DEPARTURES 
FROM ZERO OR 
EQUILIBRIUM 
GRADIENTS 


•SUDDEN APPLICATION OR 
REMOVAL OF LARGE 
GRADIENTS INWALL 
TEMPERATURE . 
PRESSURE, WALL 
INJECTION 

•TRANSITIONAL FLOWS 


•HIGH REYNOLDS NUMBER ‘LOW REYNOLDS NUMBER 
L ZEROTH ORDER METHODS 


MEAN FIELD METHODS 


MEAN TURBULENCE FIELD METHODS (WITH LENGTH SCALE EQUATION) 


In the present review an operational definition of equilibrium is used, taken 
from reference 46: [equilibrium refers^] "to a layer that has completely settled 

down after a discontinuity and is developing with its new boundary conditions 
with no 'memory' of the discontinuity." Succeeding portions of the present 
report cover these three closure approaches in considerable detail, especially 
the mean field approaches. 
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FLOW PHENOMENA PECULIAR TO CALCULATION OF COMPRESSIBLE 


TURBULENT BOUNDARY LAYERS 


Calculation procedures for compressible turbulent boundary layers are 
based, to a great extent, upon techniques, modeling constants, etc., developed 
originally for the low-speed case. This section discusses many of the signifi- 
cant differences and new or altered physics which occur in the compressible 
case, as compared with the low-speed situation. The discussion is intended to 
aid in evaluating the applicability of low-speed results to the compressible 
(particularly high Mach number) case and to indicate possible pitfalls and 
sources of inaccuracy in the calculation of compressible turbulent boundary 
layers . 


Normal Pressure Gradient 

Shown in figure 1 are typical static-pressure distributions measured across 
high Mach number turbulent boundary layers (two cases are shown, helium (ref. 47) 
and nitrogen (ref. 48)). In both cases the boundary layer involved was on a 
nozzle wall with only a small local longitudinal pressure gradient, although 
there are significant free-stream static-pressure variations due to uncanceled 
Mach waves. (See ref. 48.) Two points are obvious from this figure: the static 

pressure is not constant across these high Mach number turbulent boundary layers 
(3p/3y = 0 for M -*■ 0), and the wall pressure value is greater than the edge 
value by approximately 50 percent for these cases (M^ 20). 

A possible origin of at least a portion of this nonconstant p(y) can be 
readily seen from a simple examination of equation (3) (normal momentum equation) 


9£ « - <L(pv' 2 ) (7) 

3y 3y 


or 


p + p v’ 2 88 Constant (8) 

Evaluating this expression between the region of peak v' 2 (y/6 = 0(0.1) and 

Pp ~ p e from fig. 2) and the wall, one obtains 

Pw~ Pe + f Pe u e 2 (9) 

u e^ 

where Pp =® f p e or 

^ i + f y m 6 2 Ip! do) 

p e u e 2 

Limited data (e.g. , ref. 49) indicate that the nondimensional velocity fluctua- 
tion levels in boundary layers are not significantly influenced by Mach number, 
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and therefore a reasonable value for \Jv'^/u e ^ of 0.06 was used (from low- 
speed data). Using f n 0.2 and y H 1-5, equation (10) becomes (to correct 
order of magnitude) 

?£ z 1 + (10“3) M e 2 (11) 

Pe 

and indicates a possible dominant influence of Mach number upon the nonconstant 
p(y) phenomenon, although the influence may not be as strong as M 2 due to f 
being a function of M. Shown in figure 2 (from ref. 50) is a comparison 
between this simple prediction (eq. (11)) and most of the available data. The 
relatively good agreement between the present simple expression (eq. (11)) and 
the available data (fig. 2) may indicate that at least some of the nonconstant 
p(y) effect is caused by the turbulence field. This is probably aggravated at 
high Mach number by the fact that the dynamic pressure associated with the turbu- 
lent fluctuations becomes a significant fraction of the static pressure level at 
high Mach number. Research by Finley (ref. 51) indicates that most of the 
Pw > Pe Problem is due primarily to inviscid disturbances whose detailed influ- 
ence is modified by the turbulence effect just described. 

Calculation experience indicates that the large p changes associated with 
static temperature variations for high Mach number boundary layers (a change in 
density by a factor of approximately 100 for M ~ 20, T w ->• T^, and y = 5/3) 
greatly override the P w /p e > 1 effect and therefore, while interesting, the 
p w /p e > 1 effect is not currently considered first order for calculating high 
Mach number boundary-layer flows. 


Influence of Compressibility and Density Fluctuation Terms 

There are no definitive, detailed measurements of the complete second-order 
correlations in a highly compressible flow (including p' and p’ terms) with 
which to assess, in a straightforward manner, the influence of compressibility 
and density fluctuation terms upon closure models used in the calculation of com- 
pressible turbulent boundary layers. There is evidently a true Mach number 
effect on turbulence structure for free flows (ref. 52), at least for the free 
shear layer case at high Reynolds number with a sustained Mach number difference 
across the shear layer. However, no important compressibility (p', p') influ- 
ence has ever been isolated for the boundary-layer case (except for p'v'; see 
discussion following eq. (6)). Detailed consideration of this area is beyond 
the scope of the present paper; the purpose here is merely to warn the reader 
that the following arguments concerning the relative absence of noticeable com- 
pressibility effects are deductive. There does not yet exist a definitive set 
of measurements to completely lay to rest the question of possible compressibil- 
ity effects upon turbulence structure, although Morkovin's arguments (ref. 26, 
pp. 367-380), which were based upon hot-wire data up to M ~ 5 where p' 
effects were neglected, have proved thus far to be correct. 

Fluctuating Mach number .- One method of evaluating the possible influence 
of compressibility upon turbulence structure is to examine the magnitude of the 
fluctuating Mach number M ’ ; that is, 
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( 12 ) 


M* % ui _ Ma 1 
a a 

For the case of maximum p' (i.e., T w -*• T t>e and T t (y) = Constant) near the 
edge of the sublayer 

ul = ul M p \/5 

a u e VT 
and for y = 1 . 4 

u 1 g-, u 1 ^e 

5 u e ^/l + M e 2 

which is only of the order of 0.2 for high M e _ and u'/u e —0.1. Therefore the 
main term in equation (12) is the second one (M a'/a), which, for a'/a = 0(0.1), 
can be of the order of 1 or greater for high M e . Therefore, for the high hyper- 
sonic case the fluctuating Mach number can be of the order of 1 and compressibil- 
ity effects may become important for accurate turbulence modeling. 

Presence of p 1 terms in Reynolds stress expression .- In equation (2) the 

complete Reynolds stress term is (pv)'u'. Expanding this term (ref. 13) one 
obtains 


(13) 


(14) 


(pv) 'u' = p u'v' + v p'u' + p'u'v' (15) 

Using the usual order of magnitude arguments, v « p and p' « p, 


(pv) 'u' p u'v' (16) 

within an accuracy of approximately 20 percent. Calculation experience (e.g., 
ref. 53) indicates that the p' terms in the Reynolds stress equation (eq. (15)) 
are only important for the case where e/y is relatively small (e/y < 100). 
Therefore, for the calculation of compressible turbulent boundary layers the 
Reynolds stress term generally assumes the sam e form as the low-speed case 
(p u'v' ) . This does not imply, however, that u'v' can be modeled in the same 
manner (or has similar values) as in the low-speed case. 


Comparison of Mean Static Temperature and Turbulent Kinetic Energy 

As will now be shown, there is a possibility that, for high Mach numbers, 
the static temperature can be of the same order of magnitude as the temperature 
equivalence of the energy associated with the turbulent velocity fluctuations. 
The basic problem is most easily seen by beginning with the expression for 
instantaneous total enthalpy 

H = c D T t = c d T + + w£ (17) 

p z p 2 2 2 
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The use of Reynolds time averaging (and recognizing that w 
dimensional boundary layers) yields 


0 for two- 


n2 


2 P T t = c p T + + 


r,2 


W' 


( 18 ) 


t ^ T t - U 2 - a ' 2 + v' 2 + w' 2 + Y 2 (19) 

2c p 2c p 2Cp 2c p 2Cp 

^ v ' 

I II III 


In order to determine the static temperature from a boundary-layer solution, 
some form of equation (19) is generally employed, using terms I and II only. 

Term III has been almost universally neglected. For M e 22, T w ^t , e » anc * 
Miocal ~ 10 (i.e., at y/6 ~ 0.5) in air, the usual terms (I and II) yield 

T/T^ « 0.0476. Assuming u'/u e ~ 0.05, the additional (high M) term (III) 
yields a value of approximately 0.0075, or approximately 16 percent of the value 
of terms I and II combined. Similar comparisons over the Mach number range are 
shown in figure 3- 

What this discussion indicates is that for very high Mach numbers, 

Tt = u 2 /2c p and T s t a tic can become of the same order as the temperature equiv- 
alence of the energy associated with the turbulent motions. Since the intensity 
of the turbulence is not generally known to within an accuracy of approximately 
20 to 30 percent and T s t a tic can depend upon the square of the intensity 
(eq. (19)), the accuracy of the computed mean density field for extremely high 
Mach number turbulent boundary layers could be quite poor. 


Precursor Transition Effect 

The precursor transition effect is characterized by the existence of large- 
scale disturbances (and subsequent breakdown into turbulence) in the outer 
region of compressible boundary layers far upstream of the nominal wall transi- 
tion point. (See fig. 4 for typical schlieren photographs of this behavior.) 

This phenomenon is quite commonly observed (e.g., refs. 54 to 56). The spreading 
rates and possible structure of these disturbances were examined in reference 57 
(see fig. 5, taken from ref. 57), while an attempt to model their influence in a 
mean field closure procedure is reported in reference 58. The major impact of 
this phenomenon upon the calculation of turbulent boundary layers is that at the 
nominal wall transition location, the outer portion of the profile is already 
transitional and turbulent in nature (see fig. 6, taken from ref. 58), and there- 
fore, the usual procedure of starting the calculation at the beginning of the 
wall transition location with a laminar profile cannot be followed. 

A simple method of including this phenomenon is to start the calculation 
upstream of the wall transition location (with a laminar-like profile) where the 
transition bursts first initiate. (See ref. 59.) Figure 7 shows a possible fur- 
ther manifestation of this effect - an increase in surface heating upstream of 
the nominal wall transition point as usually defined (taken from ref. 60). As a 
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result of the outward movement of the critical layer with Mach number (fig. 8 
(from ref. 61) for near adiabatic conditions), this precursor effect occurs fur- 
ther upstream and becomes more important as Mach number increases. Unless a 
measured starting profile is available, accurate compressible-boundary-layer cal- 
culations, at least for the low to moderate Reynolds number range, should include 
some consideration of this precursor influence. (See also ref. 32.) Experimen- 
tal evidence indicates that the magnitude of the precursor transition effect may 
be influenced by T w /T e , the effect perhaps becoming less pronounced at lower 
values of T w /T e . 


Large p' Levels 

The usual assumption in second-order closure approaches is that most 
pressure fluctuation terms, especially the so-called pressure-dilitation terms 
(ref. 52) can be neglected. The reasons for this assumption are threefold: 

(1) The limited data for low speeds indicate that these terms are indeed small; 

(2) since M* is usually small, the influence of compressibility (which is pri- 
marily represented in the second-order approaches by the p' terms) must be 
small; (3) there is an absence of data concerning the actual magnitude of these 
terms. Although no data for the terms themselves exist, there is limited data 
for p' itself, which will now be examined. 

Figure 9 indicates the level of rms wall-pressure fluctuations (normalized 
by the external mean static pressure) as a function of Mach number (from 
ref. 47). These data indicate that for moderate hypersonic Mach numbers 
(M 8 to 12) the wall-pressure fluctuation levels approach 10 percent, which is 
a fluctuation level typical of the longitudinal velocity field. In fact, data 
from reference 62 at M e = 9-4 indicate that the turbulence is dominated by 
high-frequency pressure fluctuations. 

The distribution of p' with distance away from the wall is shown in fig- 
ure 10 (from ref. 62). These data indicate that p' levels are large, not only 
at the wall, but across the sensibly turbulent portion of the boundary layer. 

The p' intensity diminishes only in the outer (intermittent) portion of the 
flow. Therefore, the p' levels for hypersonic turbulent boundary layers can 
evidently be large; this calls for a close, careful examination of second-order 
closure schemes for the high M case (and also for separated flows, where p' 
levels are large e ven at low Mach number). Of especial importance are p' 
terms such as p'(3u|/3x^). (See also ref. 63-) 


Definition of Boundary-Layer Thickness 

Difficulty in defining the boundary-layer thickness occurs because of large 
differences between and 6 velocity f ' or hi 6 h M boundary layers. A 

typical result for M 20 (from ref. 47) is shown in figure 11, where the nom- 

inal <5 velocity is approximately 50 percent of Sdensity (° r ^pitot)- Th e 
pitot edge is easier to measure but there are two fundamental questions: with 

a relative absence of mean shear above 5 V gXocity> what is the turbulence shear 
stress (and turbulence eddy) content in this region, and what thickness does one 
use to scale conventional mean field closure models (such as mixing length) for 


13 


the high M ease? These questions also involve the fact, as previously dis- 
cussed, that the density becomes difficult to determine accurately for high M 
because it is a small difference between two large numbers. These questions are 
still open, but limited information (for instance, ref. 64) indicates that 
^velocity should be used as the scaling function. 


Energy Loss Via Acoustic Waves 

This subject has received little attention since the original work of 
Laufer (ref. 26, pp. 381-393). The physical problem arises because of the 
increasing intensity of turbulent boundary-layer sound radiation with increasing 
Mach number. For large enough radiated sound levels the amount of energy car- 
ried by these waves can be appreciable and represents a new energy dissipation 
mechanism in high Mach number compressible flows. 

In reference 26 Laufer computed that for M e « 5, the energy radiated away 
is approximately 1 percent of the work done by turbulent shear and therefore is 
negligible. However, if Laufer’ s equations hold for the case with M e = 20, 
T w /T t -*■ 1, and y = 5/3 , the same calculation yields radiated energy of 25 per- 
cent which is no longer negligible. For the nozzle wall boundary layers or for 
boundary layers measured on wind-tunnel models which have turbulent boundary 
layers, the flow probably reaches an equilibrium state involving a balance 
between the absorption of acoustic energy radiated by the turbulent flows which 
surround the local flow and the acoustic energy radiated by the local flow 
itself, that is, a balance between gain and loss as far as acoustic energy is 
concerned. This balance probably results in a lower net loss and thus may tend 
to obscure (in ground facilities) the true importance of this energy radiation 
effect. Only measurements in a "quiet tunnel" (ref. 65) or on a flight vehicle 
can clarify the true influence of energy loss via acoustic waves. 

There are other problems peculiar to the calculation of compressible turbu- 
lent boundary layers such as (1) persistence and importance of wall temperature 
history effects, (2) increasing predominance of low Reynolds number amplifica- 
tion, and (3) variable edge entropy, but these topics are more conveniently dis- 
cussed in connection with the mean field closure methods. It should be noted 
that the nondimen sional burst period for compressible flows is approximately the 
same as that for low-speed flows (ref. 66). 


TRADITIONAL PREDICTION METHODS 

Methods for calculating turbulent , compressible boundary-layer development 
have been a topic of extensive research for many years. Until the advent of the 
high-speed computer, most of this research effort was concerned with prediction 
methods which could be applied by using only hand calculations or very simple 
machine computations. Even now, work is continuing on these easily calculable 
methods because of their inherent accuracy (due to extensive empirical valida- 
tion) and because of the ease and speed with which these methods can be applied 
to engineering problems. Traditional methods can be classified into two general 
categories: integral methods and correlation methods. Numerous survey papers 

are available (e.g., ref. 18 and pp. 181-229 of ref. 28) which attempt to evalu- 
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ate the accuracy and limits of applicability of various traditional methods for 
compressible two-dimensional turbulent flows. 

The objectives of this section are to describe and discuss the merits of 
each category of traditional methods. In addition, several illustrative exam- 
ples of each category will be discussed and compared with experimental data. 
Comments concerning the accuracy and applicability of these prediction methods 
will be made. A discussion of compressibility transformations is first pre- 
sented since transformations are frequently used in both correlation methods and 
integral methods; next, a discussion of correlation methods including compressi- 
ble law-of-the-wall/law-of-the-wake formulations is given; finally, integral 
methods are discussed. Integral methods may use a transformation technique 
and/or a correlation method in the complete solution procedure. 


Transformations 

One of the earliest approaches to predicting compressible turbulent flows 
was to seek a transformation which, when applied to the governing equations for 
compressible flow, yields identically the incompressible equations. This tech- 
nique has been used with success for laminar boundary- layer flows which satisfy 
certain requirements on the relationship between kinematic viscosity' and tempera- 
ture. (Reynolds stress terms do not appear in the laminar boundary-layer equa- 
tions.) The obvious advantage of this approach for turbulent flows is that, if 
successful, the more extensive knowledge of the mechanism of turbulent momentum 
transfer for incompressible flows can be used to predict compressible flows. 

The transformation concept assumes that the companion incompressible flow result- 
ing from the transformation is physically observable; in addition, it is usually 
necessary to assume an invariance hypothesis concerning some transformation 
scale. Beckwith (ref. 18) points out several examples where the physical con- 
cepts of a compressibility transformation are violated because of a lack of cor- 
respondence of the transformed incompressible flow. These include the following: 
(1) Normal pressure gradients which may be important in high-speed flows do not 
exist in incompressible flow and, therefore, are not transformed; (2) dissipa- 
tion effects in heating calculations are not properly transformed; (3) large nor- 
mal temperature gradients in the transformed incompressible flow result from the 
transformation and, according to the equation of state, cannot exist in the 
constant-density flow; (4) fluctuating density terms in the compressible formula- 
tion have no counterpart in the corresponding constant-density flow. 

The earliest transformation for turbulent flow was presented by Dorodnitsyn 
(ref. 67) who considered only the Von Karman momentum integral equation. Mager 
(ref. 68) was the first to attempt to transform the partial-differential equa- 
tions for turbulent boundary-layer flow. Examples of other attempts to define a 
transformation in the same time frame were Culick and Hill (ref. 69), Burggraf 
(ref. 70), and various reference temperature or enthalpy methods for zero- 
pressure-gradient flows (e.g., refs. 71, 72, and 73)- 

Some years later, Coles (ref. 15), criticizing the assumption made in ref- 
erences 68 and 70 of invariant turbulent shear and stream function under the 
transformations, proposed a more physically acceptable transformation in which 
the adiabatic, compressible, and the constant-density flows are assumed to be 
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related by three scaling parameters a(x), ri(x), and £(x). The first param- 
eter relates the stream functions of the two flows, the second is a multiplica- 
tive factor of the Dorodnitsyn-Howarth scaling of the normal coordinate, and the 
third relates the streamwise coordinates of the two flows. An additional assump- 
tion pertaining to the invariance of a Reynolds number characterizing the law-of- 
the-wall region of the boundary layer is necessary to complete the transforma- 
tion. This assumption, which Coles called the substructure hypothesis, provides 
a substitute for a reference state utilized with many theoretical approaches. 
Coles' transformation was then modified or extended as follows: (1) Crocco 

(ref. 74) modified Coles' transformation to include effects of heat transfer 
and pressure gradient; (2) Baronti and Libby (ref. 75) replaced Coles' substruc- 
ture hypothesis with a sublayer hypothesis based on experimental observation; 

(3) Jeromin (ref. 76) extended Coles' transformation to include effects of mass 
transfer; (4) Lewis, Kubota, and Webb (ref. 77) defined new coordinates consis- 
tent with Coles' transformation and eliminated the need for a substructure or 
sublayer hypothesis (they found, however, that dissipation effects for compressi- 
ble flows invalidate the transformation for high wall temperatures); (5) Economos 
and Boccio (ref. 78) empirically modified Coles' law of the wall/law of the 
wake and introduced two compatability equations which provide the closure condi- 
tions. Coles' approach and the subsequent companion work represent the primary 
advances in transformation theory in recent years. 

Comparisons with experimental data have shown that transformations yield 
good results only for moderate Mach numbers (M g 6) and for moderate wall heat- 
ing (T w /T^ £ 0.5) with zero or mild pressure gradients (refs. 75, 77, 79, and 
80). (In ref. 79, it was observed that while compressibility had little effect 
on mixing length for flat-plate-type turbulent boundary layers, Mager's trans- 
formation predicts a large effect.) By empirically modifying Coles' transforma- 
tion, Economos and Boccio (ref. 78) were able to extend the range of agreement 
with experimental data beyond any other method. While this empirical modifica- 
tion was not physically appealing (transformation of impermeable wall case 
yields mass transfer in transformed plane), agreement with experiment was sub- 
stantially improved. This is illustrated in figure 12 where the methods of 
Baronti and Libby (ref. 75) and Economos and Boccio are compared with experi- 
ment; here, local skin friction values were obtained from experimental velocity 
profiles according to each transformation approach (by fitting law-of-the-wall 
profiles to experimental profiles) and normalized by a reference value. (See 
ref. 80 for details.) This reference value is equivalent to the measured skin 
friction as illustrated in figure 12. As the ratio of wall temperature to total 
temperature is decreased, a systematic error appears in the skin friction 
obtained using Coles' transformation with the sublayer hypothesis as recommended 
by Baronti and Libby. The empirical modification of Economos and Boccio appar- 
ently eliminates this disagreement. 

The ability of Coles' original transformation to predict flat plate skin 
friction and heating for Mach 4 to 13 and 0.1 4 1 T w /T^ ^0.7 is illustrated in 
figure 13- The experimental data are transformed according to Coles' method and 
compared with a good incompressible prediction (ref. 81). The prediction of 
skin friction is generally poor; this poor agreement results from an effect of 
T w /T^ not accounted for in the transformation. (See ref. 82.) The improved 
heating prediction in figure 13 reflects the use of a Reynolds analogy factor 
( 2 N 5 |VCf) to adjust the level of the data. 
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Both references 77 and 78 applied their versions of Coles' transformation 
with an incompressible integral method, and examples of the results are shown in 
figure 14. Here, the skin friction and momentum thickness for the waisted body 
of revolution from reference 83 are compared with predicted values from each 
method. These comparisons properly belong in a later section on integral meth- 
ods but are shown here to illustrate the accuracy of these transformation meth- 
ods. The agreement of the method of reference 78 (fig. 14(a)) with experiment 
is good and rivals some of the more advanced numerical methods, while the predic- 
tions of reference 77 (fig. 14(b)) are obviously poor and unacceptable for engi- 
neering purposes. (However, curvature and laminar ization effects may be present 
in some of the data of ref. 83 and thus could influence these comparisons.) 

These comparisons of the results from transformation methods with experi- 
mental data and other comparisons in the cited references lead to the general 
conclusion that while some success for adiabatic flows is evident, transforma- 
tion methods are not presently desirable for general applications, a conclusion 
also reached in references 18 and 28 (pp. 181-229). 


Correlation Methods 

Correlation methods are perhaps the most popular of the traditional methods 
because they are generally simple to apply and some are quite accurate. Even 
though these methods are strictly applicable only for two-dimensional or axisym- 
metric flows with zero pressure and wall temperature gradients, they are widely 
used for parametric design studies and preliminary estimates. The empirical 
input to some of these methods enhances the accuracy and allows easy modifica- 
tion of the method to include new empiricism. Some of these methods provide pre- 
dictions for wall shear stress only, while others also allow specification of 
the boundary-layer profiles. Additional assumptions for the Reynolds analogy 
factor 2Ngt/Cf and the recovery factor are required in order to calculate wall 
heat transfer for those methods which only predict wall shear. 

Spalding and Chi (ref. 81) have presented an excellent summary of available 
correlation methods for compressible turbulent flow up to 1964. Since 1964 the 
only correlation method to arise is that of White and Christoph (ref. 84). 
According to Spalding and Chi, correlation methods may be classified as follows: 

(1) Methods using Prandtl or Von Karman differential equations (i.e., 

Prandtl or Von Karman mixing length concepts) 

(2) Theories based upon other differential equations 

(3) Theories based upon a fixed velocity profile 

(4) Theories based upon incompressible formulae with reference properties 

In categories 1 and 2 the shear stress is assumed to be constant through 
the boundary layer and equal to its wall value; in category 3 the velocity pro- 
file is assumed indpendent of compressibility; in category 4 the incompressible 
relations are assumed to apply for compressible flows if the gas properties are 
evaluated at a reference temperature or enthalpy where the reference temperature 
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is a function of Mach number, ratio of wall to edge temperature, and recovery 
factor. Some examples of the more popular methods in each category are as 
follows: 

Category 1 - Van Driest I (ref. 13) and Kutateladze and Leont'ev (ref. 85), 
which use a Prandtl mixing length, and Wilson (ref. 86), Van Driest II (ref. 87), 
Harkness (ref. 88), Deissler and Loeffler (ref. 39), and Moore (ref. 89), which 
use a Von Karman mixing length. 

Category 2 - Li and Nagamatsu (ref. 90) and Kosterin and Kashmarov 
(ref. 91). 

Category 3 - Cope (ref. 92) and Monaghan (ref. 93). 

Category 4 - Summer and Short (ref. 73) and Eckert (ref. 72). 

In addition to these methods, Spalding and Chi (ref. 81) present a method which 
uses a Van Driest type analysis to determine a function F c and empirical data 
to determine a function Fp. In their method the assumption is made that F c Cf 
varies uniquely with FpR according to an incompressible law (R is a Reynolds 
number) and that F c and Fp are functions only of Mach number, ratio of wall 
temperature to edge temperature, and recovery factor. 

Numerous survey papers (refs. 16, 80 to 82, 84, and 94 to 98) have 
attempted to determine which’ of the correlation techniques best predicts avail- 
able experimental data for compressible turbulent boundary-layer flow. The 
method which is found to be the "best" in each of these survey papers is usually 
either Van Driest II (ref. 87), Spalding and Chi (ref. 81), or a reference tem- 
perature approach (e.g., Eckert, ref. 72). The choice of a "best" method is 
influenced in each of the survey papers by particular selections of virtual ori- 
gin, Reynolds analogy factor for heating calculations, and recovery factor. 

(See ref. 82.) 

Comparisons of experimental skin friction and heat transfer data with pre- 
dicted values obtained from the methods of Van Driest (ref. 87), Spalding and 
Chi (ref. 8l ) , and Eckert (ref. 72) for Mach numbers 4 to 10 and ratios of wall 
temperature to total temperature from 0.14 to 0.7 are shown in figures 15 to 17. 
The data were obtained in wind tunnels on flat-plate models with dp/dx = 0 and 
dT w /dx = 0; skin friction was measured with balances and wall heating was mea- 
sured by transient techniques. The comparisons are presented in the form F c Cf 
against Fp Rg or Fp R x vo where F c Cf is equivalent to the incompressible 
value of skin friction x an<5 Fp Rq or Fp R x vo is equivalent to the correspond- 
ing incompressible value of Reynolds number. (See ref. 82 for further explana- 
tion . ) In this manner , all the transformed experimental data can be compared 
directly with values from a good incompressible skin friction law (law from 
ref. 81 used herein) to judge the efficacy of each method. 

The comparison of the data with predicted values from the Van Driest II 
method (ref. 87) is shown in figure 15. For this comparison, a momentum thick- 
ness Reynolds number is used, and a Reynolds analogy factor of 1.0 is assumed, 
as recommended for best results with the Van Driest II method by Hopkins and 
Inouye in reference 98. The comparison is favorable for heat transfer, but the 
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transformed skin friction data fall generally below the prediction. As reported 
in reference 82, the prediction error also appears to be a function of VTf 
A similar comparison for the method of Spalding and Chi (ref. 81) is shown in 
figure 16. Here, the virtual origin location was assumed to be near the end of 
boundary-layer transition, and Von Karman's Reynolds analogy factor (used as sug- 
gested in the appendix of ref. 95) was applied as recommended for the Spalding 
and Chi method in reference 82. The prediction of both the heating and skin 
friction data in figure 16 is quite good. The results from Eckert's method 
(ref. 71) are shown in figure 17 using the same virtual origin and Reynolds anal- 
ogy factor as previously used with the Spalding and Chi method. The comparison 
between data and prediction is again very good in figure 17, but there appears 
to be a discernible variation of prediction errors with T w /Tt- (See ref. 82.) 
Considering the wide range of flow conditions of the data, each of the three 
approaches provides a credible prediction method for zero-pressure-gradient 
turbulent-boundary- layer flows. 

Comparisons of predictions from Eckert's method (ref. 72) with data from 
flight are shown in figures 18 and 19. Data obtained on conventional aircraft 
such as the A-5A, the Mirage IV, and the XB-70-1 (refs. 99 to 101) are shown in 
figure 18. The agreement with Eckert's method is good in each case. Heating 
data from sharp and blunt cones in flight are shown in figure 19 (from ref. 96) 
where Colburn’s Reynolds analogy (ref. 102) was used for the calculation and the 
virtual origin was assumed to occur at the stagnation point or cone tip. A 
cone/flat plate transformation after Van Driest (ref. 103) was applied. The 
agreement between prediction and data is surprisingly good over the wide range of 
data included (3 g M g 13, 0.2 g T w /T e g 2.3). It appears therefore that these 

correlation methods provide accurate predictions of turbulent skin friction and 
heat transfer for a wide range of flow conditions (for dp/dx ~ 0, dT w /dx ~ 0) 
in flight as well as in wind tunnels. 

Compressible law of the wall .- The law-of-the-wall velocity profile has 
been useful in arriving at quantitative analytical results for incompressible 
flows. This incompressible law is empirical in nature but can be derived 
directly from mixing length concepts; that is, 

T ■ 5l2 (l jf 

and 

l = Ky (Prandtl) 


so 



| Integration yields (p = Constant and T = x w = Constant) 

I L = u + = 1 In y + + C 

u T K 


( 20 ) 


( 21 ) 
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where (tt/K)w is a wake function which accounts for the wake-like behavior of 
the outer portion of the turbulent boundary layer. Alternately, by evaluating 
equation (22) at the edge of the boundary layer and subtracting the resulting 
equation from equation (22), the velocity defect equation is obtained 

Ue ~ U = - 1 In 2 + 5(2 - w) (23) 

u T K 6 K 

An expression for the incompressible profile from the wall to the outer 
edge of the buffer layer (beginning of fully turbulent region) was presented by 
Spalding (ref. 105) and Kleinstein (ref. 106). Modifications to the law of the 
wall for surface roughness (e.g., ref. 107) and wall injection (e.g., ref. 108) 
are available. Equation (23), the velocity defect distribution, is independent 
of surface roughness. The law-of-the-wall/law-of-the-wake expressions have been 
used extensively and with great success in studies of turbulent boundary layers 
either to obtain such parameters as Cf and <5 from experimental velocity pro- 
files or as an auxiliary equation in integral or numerical prediction methods. 

The successful application of the law-of-the-wall/law-of-the-wake velocity 
profile for incompressible flows has naturally led to numerous attempts to apply 
similar concepts to compressible flows. A brief description of some of the more 
prominent of these studies follows. 

Van Driest (ref. 13) in 1951 derived a compressible law of the wall for 
zero-pressure-gradient flows by solving the turbulent-boundary-layer equations 
with Prandtl's mixing length formulation and a laminar and turbulent Prandtl num- 
ber of unity. By comparing the resulting law of the wall from Van Driest' s anal- 
ysis with the incompressible law of the wall, it is clear that the effects of 
compressibiliy can be accounted for by defining a generalized velocity u* in 
place of u. (See Maise and McDonald (ref. 79).) This generalized velocity is 
defined as 

u* = u 1 sin~1 2A^(u/u rr> ) - B _ f /P du 

001 L(B 2 + A 2) 1/2 J J IpwJ 

where p is from Crocco ? s relation and 

t? = [(y - n/alMj 
T »/T„ 
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B = 1 + C(Y - D/2]M J . , 

More recently Maise and McDonald (ref. 79) applied these concepts by assum- 
ing that Coles' law of the wall was valid for compressible flow if the veloci- 
ties are defined as the generalized velocity. The resulting equation for the 
velocity defect was 

U I ~ u * = - 1 In y + C(2 - w) 
u x K Z 

where w is Coles' tabulated wake function (ref. 104). Using Clauser's results 
(ref. 1) for incompressible flow to define the constants, the final expression 
is 

* 

u ~ - u * = -2.5 In 1 + 1.25(2 - w) (24) 

u t 6 

Maise and McDonald found that experimental velocity profiles for adiabatic, com- 
pressible turbulent boundary layers from Mach 1.5 to 5 through a wide range of 
Reynolds numbers were well correlated by equation (24). However, the velocity 
profiles for nonadiabatic flows in the same range of Mach number were poorly 
correlated. 


Fernholz (ref. 109) used a similar analysis with generalized velocities but 
did not use Clauser's constants. Instead, experimental data were used to corre- 
late the constant as a function of the ratio of wall temperature to adiabatic 
wall temperature and Reynolds number. The formulation was 


H! = 1 in + Fl (25) 

U'j* K 

where 

p, _ f/^aw " T„ p e u e 9^ 

1 • TT5— ’ — j 


Fernholz also included the velocity defect law but redefined the 6-coordinate 
as was done for incompressible flows. The resulting equation was 


# 

u - u* 

QQ 

Ut 


-M In L - N 

A* 


( 26 ) 


where 



A A 

M = 4.7 and N = 6.8 for adiabatic walls. Equation (25) provided a good corre- 
lation of experimental velocity profiles for Mach 5 to 8 with moderate heat 
transfer. Equation (26) also correlated the defect region of these profiles 
for Rq = p e u e 9/y w greater than 1000 to 1500. 
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Mathews, Childs, and Paynter (ref. 110) using generalized velocity concepts 
fitted the velocity defect equation (eq. (24)) to experimental velocity profiles 
(by the method of least squares) and thereby obtained the skin friction and 
boundary-layer thickness which allowed the best profile fit. The test cases 
involved both normal and conical shock interactions with turbulent adiabatic 
boundary layers at supersonic speeds as well as a flat-plate supersonic flow. 

The velocity profile fits were good, and the resulting skin friction and 
boundary-layer thickness values compared favorably with experiment and other 
predictions . 

Sun and Childs (ref. Ill) modified the analysis of reference 110 by using a 
more realistic shear stress distribution through the boundary layer. Instead of 
assuming T □ T w = Constant, the expression 

k - ' - (if <27> 

was used. Equation (27) with b = 1 is a reasonable fit to experimental data 
for both subsonic and supersonic flows (ref. 25). This modification satisfies 
the physical condition that the velocity gradient be zero at y = <5 and pro- 
vides more accurate values of boundary-layer thickness. The profile fit and the 
resulting skin friction were little changed from the analysis of reference 109. 

Kane (refs. 112 and 113) used the compressible counterpart of equation (1) 
which is 

u ++ = 1 In y + + C (28) 

K 

where 

/2 du+ (29) 

Equation (29) was evaluated by using the Crocco temperature profile and Squire's 
definition of t/t w (ref. 114) which includes the effects of wall injection. 
Inclusion of an empirical extension of Kleinstein's buffer layer profile 
(ref. 106) to compressible flow as well as Coles' wake parameter completed the 
formulation of the compressible law of the wall/law of the wake. Kane finds 
that three profile functions remain undefined in his compressible law of the 
wall/law of the wake which must be obtained from a fit of experimental velocity 
profile data. A multiple regression analysis is then used to define the best 
functional variation of these three functions with appropriate fluid variables. 

Squire (ref. 114) used an approach similar to Kane's method to extend the 
compressible law of the wall to flows with wall injection. By fitting experimen- 
tal velocity profile data, the integration constant in the law of the wall was 
determined as a function of v w /u T and Mach number. The slope of the log 
region K was found to be independent of mass injection. Chen (ref. 115) 
extended the generalized velocity concepts used by Maise and McDonald (ref. 79) 
to flows with rough wall, heat transfer, and pressure gradient. White and 
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Christoph (ref. 1 1 6 ) included effects of longitudinal pressure gradient on the 
shear stress distribution and then, using a method similar to Squire (ref. 114) 
with v w = 0, obtained a closed form solution for the law-of-the-wall velocity 
profile. 

Danberg (ref. 115) used Coles' law-of-the-wall/law-of-the-wake profile and 
a reduced velocity 



to account for compressibility effects. Four profile parameters were determined 
from least square fits of Danberg 's compressible law of the wall to experimental 
velocity profile data from Mach 2 to 6 for adiabatic wall conditions. Danberg 
attempted a similar definition of a law of the wall for the temperature profile, 
but a paucity of experimental data prevented definitive results. 

This discussion of the status and development of a compressible law of the 
wall indicates that reliable expressions for a compressible law of the wall 
exist for adiabatic wall conditions and Mach numbers less than 6. However, for 
flows with significant heat transfer, the definition of a compressible law of 
the wall is still unclear. The unresolved question of the effect of heat trans- 
fer on the applicability of a compressible law of the wall is clearly illus- 
trated when the results of Maise and McDonald (ref. 79), which show poor corre- 
lation of compressible velocity profiles with heat transfer, are compared with 
the results of Gran, Lewis, and Kubota (ref. 118) which show good correlation of 
compressible velocity profiles with pressure gradient and heat transfer, both 
studies using nearly the same compressible law-of-the-wall formulation; the dif- 
ference between the two formulations was that A* (see eq. (26)) was used 
instead of 6 in the velocity defect law in reference 117. It thus appears 
that further progress in defining a general law-of-the-wall velocity profile for 
compressible turbulent boundary layers will be paced by progress in obtaining 
detailed and accurate experimental data through a wide range of flow variables. 


Integral Methods 

All integral methods solve the (Von Karman) integral momentum equation 
along with various auxiliary relations. The two-dimensional Von Karman momentum 
integral equation is obtained by integrating the x-momentum equation. When nor- 
mal stress terms are neglected, the following equation is obtained: 

® . L ^/h . 2 . is *•) + _J 3-(p e « - f p dy) = £f ( 30 ) 

dx u e dx ^ p e du e ) PgUe 2 dx \ J 0 ) 2 

where 



23 



If p is independent of y and the free-stream flow is adiabatic, then 

u e d Pe -m2 
— - 1 

Pe du e 

and equation (30) becomes 

d6 + 9_ ^e(H + 2 - M e 2 ) = 
dx u e dx 2 

There are two basic types of integral methods - those which solve the equa- 
tions in the physical plane and those which use a compressibility transformation 
and then solve the equations in the incompressible plane. Most of the earlier 
methods used a compressibility transformation, while more recent methods favor 
solving the equations in the physical plane. Examples of auxiliary relations 
used in the numerous methods available in the literature are as follows: 

(1) Moments of the integral momentum equation, which are obtained by multi- 

plying the momentum equation by y or u and integrating across the 
boundary layer (these moment equations contain unknowns involving tur- 
bulent shear stress or dissipation integrals which must be defined in 
terms of known quantities) 

(2) Entrainment equation and/or lag equation 

(3) Specification of velocity, temperature, and/or shear stress profiles 

(4) Specification of form factor, wall shear, and/or wall heating (Reynolds 

analogy) 

By using selected auxiliary relations, the problem is reduced to solving a set 
of quasilinear, coupled ordinary differential equations by an available solution 
technique . 

Integral methods have advantages which have assured their continued use 
through the years. Primary advantages are the following: (1) Solution proce- 

dures are fast and easily programed; (2) starting procedure is simple; (3) less 
detailed information on turbulence is necessary; (4) integration process eases 
restriction on accuracy of profile shapes. The disadvantages of such methods 
are that considerable empiricism is necessary to close the equation set (i.e., 
relying upon empirical input restricts accuracy and range of application) and 
that nonequilibrium effects are difficult to include. The latter restriction 
becomes particularly important for high Mach number flows. Of course, if a com- 
pressiblity transformation is used, the integral method is subject to all the 
uncertainties inherent in the transformation. 


24 


Most earlier integral methods utilized a compressibility transformation in 
the solution procedure. Examples of these are the methods Reshotko and Tucker 
(ref. 119), Englert (ref. 120), Mager (ref. 68), Culick and Hill (ref. 69), and 
Spence (ref. 121). More contemporary methods using transformations are those of 
Sasman and Cresci (ref. 122), Camarata and McDonald (ref. 123), Flaherty 
(ref. 124), Zwarts (ref. 125), and Kiister (ref. 30, pp. 19-1 - 19-11); exten- 
sions of Head’s entrainment method (ref. 126) to compressible flows by transfor- 
mations were given by Standen (ref. 127), So (ref. 128), and Green (ref. 129). 

In recent years, integral methods which do not use a compressiblity trans- 
formation but rather use auxiliary relations which are valid in the compressible 
plane have found considerable favor. A representative list of such methods 
includes those of Miller (ref. 130), Michel, Quemard, and Cousteix (ref. 131), 
Johnson and Boney (ref. 132), Pinckney (ref. 133), Reeves (ref. 30, pp. 6-1 - 
6-A2-2), and White and Christoph (ref. 84). In addition, Green (refs. 129 and 
134) and Green, Weeks, and Brooman (ref. 135) extended Head's entrainment method 
(ref. 126) to compressible flows also using physical variables. Obviously, dis- 
cussion of all these methods is not feasible in the present paper; therefore, a 
representative example of each type of method is discussed. The particular meth- 
ods chosen were a method using a transformation (Flaherty, ref. 124), a method 
using entrainment concepts in physical variables (Green, Weeks, and Brooman, 
ref. 135), a method specially configured for adverse gradients (generated by con- 
cave surfaces) also in physical variables (Pinckney, ref. 133), and a method in 
physical variables but using a compressibility transformation for a portion of 
the solution (Reeves, ref. 30, pp. 6-1 - 6-A2-2). The structural highlights of 
each of these four calculation procedures will be discussed along with compari- 
sons with typical experimental data. 

Method of Flaherty .- This method (ref. 124) is basically a modification of 
the procedure of Reshotko and Tucker (ref. 119) and solves the momentum integral 
and moment-of-momentum equations in the transformed plane defined by using 
Stewartson's transformation (ref. 136). The skin friction coefficient is 
obtained from the Ludwieg and Tillman incompressible equation (ref. 137) using 
reference temperature concepts (ref. 72). An empirical expression is used by 
Flaherty for the shear stress integral through the boundary layer and is the 
main improvement over the earlier Reshotko and Tucker method which used a con- 
stant shear stress. A provision for calculating wall heat transfer based upon 
the energy deficit in the boundary layer is also included. 

Comparisons of predictions from Flaherty's method with two sets of experi- 
mental data are shown in figure 20. These example data were chosen to illus- 
trate the ability of the method to predict boundary-layer growth along flat or 
curved surfaces with favorable pressure gradient. Good prediction of the 
boundary-layer thickness data from reference 138 at Mach 1.5 (fig. 20(a)) is 
achieved over the entire length of the test plate. Reasonable agreement with 
the data from reference 83 at Mach 2 (fig. 20(b)) on the waisted body of revolu- 
tion is also obtained, but the agreement of prediction with both momentum thick- 
ness and skin friction data deteriorates in the region of adverse pressure 
gradient (x/L > 0.7). Comparisons of this method with other available data 
in references 30 (pp. 181-229) and 124 indicate generally good predictions of 
boundary-layer integral properties up to Mach 6 with moderately cooled walls and 
moderate pressure gradient. 
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Method of Green, Weeks, and Brooman .- This integral method (ref. 135) is 
essentially an extension of Head's method (ref. 126) to compressible flow and 
attempts to account for the influence of upstream flow history on the turbulent 
stresses. The procedure consists of solving in the physical plane the integral 
momentum equation, the entrainment equation, and a rate equation for the entrain- 
ment coefficient. This rate equation or lag equation is developed from an equa- 
tion for shear stress derived from the turbulent kinetic energy equation. In 
addition, an empirical factor is introduced to account for the variation of inte- 
gral parameters with Mach number in equilibrium flows. The method as formulated 
is restricted to adiabatic flows. 

Two comparisons of predictions from this method with experimental data are 
shown in figure 21 . The first is a comparison with the data of reference 83 
(fig. 21(a)). The prediction is in good agreement with the experimental surface 
shear and momentum thickness, even at the rear (adverse pressure gradient) of 
the waisted body. Good agreement is also obtained with the data of Lewis, 

Gran, and Kubota (ref. 139), as shown in figure 21(b). Here the pressure gradi- 
ent is zero up to x~ 35.6 cm, adverse from x « 35.6 cm to 45.7 cm, and then 
favorable thereafter. The good agreement seen in figure 21 is also evident in 
other comparisons shown in reference 135 for adiabatic conditions up to Mach 4. 

Method of Pinckney .- This method (ref. 133) was especially tailored to pre- 
dict turbulent boundary-layer development on compression (concave) surfaces com- 
mon to hypersonic air breathing engines. The set of equations which is solved 
consists of the momentum integral equation, the moment-of-momentum equation, and 
the integral energy equation. Auxiliary relations consist of (1) assuming a 
Crocco-like temperature profile which satisfies the total energy deficit across 
the boundary layer as determined from the net heating along the surface (history 
effect), (2) assuming an empirical shear stress distribution across the boundary 
layer, and (3) assuming the Spalding and Chi (ref. 81) skin friction prediction 
(suitably shifted) is valid in pressure gradient flows (as well as the resulting 
heat transfer coefficient obtained from Spalding and Chi skin friction with Von 
Karman's Reynolds analogy factor (ref. 140)). Provisions are made to allow cal- 
culations to proceed through transition and into fully turbulent boundary-layer 
flow. 


Comparisons of Pinckney's method with experimental data from reference 1 4 1 
are shown in figure 22. These integral boundary-layer-thickness data were 
obtained on an axisymmetric compression surface for both Mach 5 adiabatic flow 
and Mach 8 cold-wall conditions. Predictions from the integral method are gener- 
ally in good agreement with the integral thicknesses for both experimental condi- 
tions shown. Comparisons with other data in reference 133 are also favorable 
and confirm reasonable accuracy for this prediction method at least up to Mach 8 
with moderate wall cooling. 

Method of Reeves .- This method (ref. 30, pp. 6-1 - 6-A2-2) solves the com- 
pressible turbulent boundary- layer equations including mass transfer by assuming 
a two- layer boundary- layer model. Mixing length concepts are assumed to apply 
in the inner layer, and this results in a compressible law-of-the-wall expres- 
sion; this expression is inserted into the boundary- layer conservation equations 
which are integrated away from the wall to a matching location. The outer (or 
wake) layer solution uses an integral momentum method, the results from which 
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are matched to the inner layer solution at a prescribed match point. The outer 
layer integral solution employs a compressibility transformation. An empirical 
expression for the shear stress integral for the outer layer is specified. 
Apparently the dynamics of the boundary-layer solutions from this method depend 
on coupling and interaction between the inner and outer solutions. 

Two comparisons of predictions from Reeves' method with experimental data 
are shown in figure 23. The predictions of the integral thicknesses for the 
data of reference 138 (fig. 23(a)) at Mach 1.5 are excellent. While the predic- 
tion of momentum thickness for the data of reference 83 (fig. 23(b)) at Mach 2.4 
is good, the skin friction data are not well predicted. Reeves found good agree- 
ment between his predictions and experimental data up to Mach 10 for a wide 
range of wall temperature ratios with and without wall mass transfer. Some com- 
putational stability problems occurred for negative pressure gradients. The 
method was computationally quite rapid. 

In summary, this brief review of various available integral solution proce- 
dures for two-dimensional compressible turbulent boundary- layer flow indicates 
that reasonably accurate predictions are possible using integral techniques for 
a wide range of flow conditions. Caution must be exercised in applying integral 
methods to flows with severe pressure gradients as well as other nonequilibrium 
effects and to flows where the empirical correlations used in the method do not 
apply. While integral methods provide fast and inexpensive calculations neces- 
sary for design and analysis of fluid systems, the limitations of any method 
must be clearly understood to prevent erroneous conclusions. As discussed ear- 
lier, transformation theory has not, thus far, developed into a tool for general 
application, and therefore, integral methods using a "complete" compressibility 
transformation should probably be avoided at present. 


NUMERICAL SOLUTION PROCEDURES 

This section considers the various alternative numerical procedures cur- 
rently being used to solve the nonlinear partial-differential equations describ- 
ing compressible turbulent boundary layers (which are parabolic in the marching 
or longitudinal direction) . There are several numerical difficulties peculiar 
to the calculation of turbulent (as opposed to laminar) compressible boundary 
layers. These problem areas include (1) presence of a thin sublayer, which 
requires either a separate wall treatment or variable grid (or coordinate trans- 
formation), (2) rapid growth of the boundary layer with longitudinal distance, 
which requires a transformation or streamline mapping procedure, and (3) alge- 
braic terms in the turbulence modeling expression which can alter the stability 
of the numerical calculation procedure. 


Solution Techniques 

The numerical solution procedures used to solve the compressible turbulent 
boundary-layer problem can be conveniently categorized in the following manner: 
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Conventional Finite-Difference Methods 

Implicit (refs. 45 (pp. 375-383), 53, and 142 to 157) 

Explicit (ref. 158) 

Finite-Difference Variants 

Shooting method (refs. 28 (pp. 551-554), 40, 159, and 1 60 ) 

Box method (ref. 1 6 1 ) 

Micro-integral (ref. 9) 

Methods Employing Analytical Functions 
Wortman approach (ref. 162) 

BLIMP (ref. 1 63) 

Method of weighted residuals or method of integral relations (refs. 164 
and 1 65 ) 

Finite element (ref. 166) 

Method of Characteristics (ref. 36) 

The conventional finite-difference methods have been the most popular. 

They usually employ Crank-Nicolson differencing and the Thomas algorithm for 
inversion of the coefficient matrix. The advantages of these procedures include 
ease of coding, numerical stability, and overall simplicity. Their disadvan- 
tages are mainly due to their essentially "brute force" approach, in that 
100 nodes in the vertical direction are usually needed to adequately represent 
a turbulent profile (even with nonuniform mesh spacing). Therefore, the proce- 
dures generally require relatively large storage and long machine time, espe- 
cially when chemical reaction effects are included. 

The shooting procedure uses finite differences in the longitudinal or 
x-direction but solves an ordinary differential equation (two-point boundary- 
value problem) in the transverse or y-direction. Several users of this approach 
have encountered serious stability problems, especially with wall heating or 
cooling and pressure gradient cases. The solutions seem to be sensitive to 
guesses of the inner boundary conditions. Only a few codes presently use the 
shooting method, and it is probably not one of the best procedures available. 

The box method is quite efficient for boundary layers and is of more recent 
vintage than the conventional finite-difference and shooting procedures. By a 
change of dependent variables the equations are reduced to a first-order nonlin- 
ear system, which is solved by Newton iteration and two-point differencing. 

This procedure has several advantages, which include high order spatial accu- 
racy, even with a rapidly varying nonuniform grid, and the small number of nodes 
required for solution (obtained by Richardson extrapolation). As a result of 
these advantages, this procedure produces solutions in 7 times less machine time 
and with much less storage than conventional finite-difference methods. 

In the micro-integral method, developed at Imperial College of Science and 
Technology, the boundary-layer growth is accounted for by use of a stream func- 
tion as a transverse variable, thus reducing the number of nodes required. Most 
versions of this procedure also use a Couette flow analysis near the wall, which 
significantly reduces the required calculation time, as the use of small nodal 
spacing near the wall to resolve the sublayer is no longer necessary. Also, in 
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this method the convective terms are integrated between grid nodes before differ- 
entiation. This solution procedure is one of the more efficient ones available, 
but the Couette flow analysis near the wall must change with different boundary 
conditions, and, since the usual formula for grid control (stream- function bound- 
ary condition) is explicit, there may be problems for large longitudinal 
increments. 

The Wortman approach is fairly recent and bears some faint resemblence to 
the shooting method in that finite differences are taken in the longitudinal 
direction only. An integrating factor is used to reduce the order of the equa- 
tions and the procedure is referred to as an iterative operator method. The 
method appears to be quite fast and has a high order of spatial accuracy. 

BLIMP (boundary-layer integral matrix procedure) has undergone extensive 
development, especially for application to flow chemistry problems. The proce- 
dure is fairly involved, and the method was specifically developed to minimize 
the number of grid nodes (especially important for the large equation set which 
results when chemical effects are included). The method is characterized by 
finite differences in the longitudinal direction, and strip integrals in the 
transverse direction with cubic spline fits and Newton-Raphson iteration. 

The method of weighted residuals (or method of integral relations) is some- 
what similar to strip integral procedures. The method uses weighting and approx- 
imating functions which must be assumed. The advantage to the procedure is 
again the small number of nodes involved, but the coding is fairly complex and 
considerable insight is often needed to select reasonable functional forms. For 
turbulent flows the computation time can be of the same order as for the conven- 
tional finite difference and the method (and the several variants thereof; 
including finite element which is a global as well as a local method of weighted 
residuals) has not been very popular. 

In the method of characteristics a solution is obtained in the outer region 
of the boundary layer only and matched to a law of the wall. For this procedure 
the viscous terms in the boundary- layer equations must be neglected. 


Systems of Independent Variables 

The simplest set of independent variables to use would be the actual nondi- 
mensional physical quantities x,y (ref. 167). An advantage of this variable 
system, when compared with the more usual transformations, is the decreased 
labor involved when changes are made in the physical specification of the prob- 
lem (i.e., one does not have to keep untransforming and retransforming for quan- 
tities or boundary conditions which are specified functions of physical vari- 
ables). Therefore, these variables are particularly useful for inclusion of 
alternate turbulence models and to determine the influence of a specified varia- 
tion of pressure in the transverse direction. However, this system of variables 
does not account for boundary- layer growth with longitudinal distance without 
periodic nodal point redistribution and is therefore seldom used in production 
codes which are expected to apply over changes in Reynolds number of several 
orders of magnitude. 
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One method of solving this boundary-layer growth problem is to use a normal- 
izing factor on y which is a function of x, such as 6, 6*, etc. This will 

keep the boundary layer within the computational grid, but longitudinal deriva- 
tives of the normalizing factors are required; if 6* is used, this quantity 
can become negative for cold walls. 

More generally, solutions are obtained using some variant of the Levy-Lees 
transformation (e.g., ref. 11) 

5 “ / Pe u e dx (31 ) 





(32) 


originally developed for laminar boundary layers. These variables reduce the 
boundary-layer growth and the influence of variable density in the computational 
domain. As a rule of thumb, n ~ 0.5 for laminar flows and n = 0.8 for tur- 
bulent flows. However, in general n = f(x) and the boundary layer may still 
grow out of (or into) the mesh (i.e., He t Constant for turbulent flows with- 
out solving for h(x)). 


Another alternate set of independent variables involves the use of a stream 
function \jj as the normal or .transverse coordinates (ref. 9). This approach 
automatically tracks the boundary-layer growth, but the boundary conditions must 
be specified or determined. One could also use the Crocco or velocity variables 
(ref. 1 68 ) where the normal coordinate is a function of u/u e . This allows an 
accurate solution with a nearly constant Au/u e step size, which corresponds to 
a highly nonuniform Ay; that is, the number of necessary nodes is very effec- 
tively minimized. However, this variable set is no t suitable for cases where 
u/u e > 1 (as transformation involves \J\ - (u/u e ) ) , and therefore, wall jet 
flows would be difficult to compute with the code. 

Two examples are given here of the possible problems and inaccuracies which 
can be encountered in numerical solutions. Figure 24 (taken from ref. 1 69 ) indi- 
cates the error in Cf associated with a change in nodal point spacing. The 
usual assumption made in numerical analysis is that the solution becomes more 
accurate as the integration interval is reduced. However, as the figure shows, 
for single-precision IBM machines with approximately 7 decimal place accuracy, 
decreasing the step size can actually increase the discrepancy between a theo- 
retical and numerical solution (as a result of roundoff error). When double- 
precision arithmetic is used (approximately 15 decimal place accuracy), the 
expected trend is obtained. The other example is given in figure 25 (from 
ref. 170). Here K = AOn+l/^Hn and is the conventional means of including var- 
iable nodal spacing in a finite-difference procedure (e.g., ref. 53). These 
results from reference 170 show that, for K i 1.0, there is some error involved 
in using variable grid spacing (at^least as used in ref. 170) and that this 
error can become appreciable for K = 1.05 (usual values of K used in solu- 
tion procedures are 1.02 to 1.04). These examples are only given to indicate 
that certain simple accuracy checks should be made on any numerical code before 
the numbers can be fully believed. 


30 


MEAN FIELD CLOSURE 


This is probably the most widely used recent closure approach for computing 
com pres sible turbulent boundary layers. In this procedure the Reynolds stress 
(p u’ v* ) is related directly to the mean velocity and density fields. This 
assumption is nearly exact for equilibrium and near equilibrium boundary layers 
where turbulence production is approximately equal to turbulence dissipation 
(e.g., pp. 275-299 of ref. 45). That is, 


y* y' 9 u cc 

3y i 

Using the Prandtl model (ref. 45, pp. 275-299) 
u»v' « We x iM 

3y 

Solving equation (34) for e and inserting this into equation (33), 
u t v « 9u oc (u'v* )3 

3 y i 4 (3u/3y)3 


(33) 


(34) 

one obtains 

(35) 


or 


CFv ' « l 2/3u\ 2 ( 3 6) 

\3y / 


which is the usual mixing-length model; that is, u'v’ is a function of mean 
velocity profile. The popularity of the mean field closure procedures is due 
to several reasons: (1) Many compressible turbulent-boundary-layer flows are 

either equilibrium or near equilibrium; (2) a wide range of boundary conditions 
can be easily and accurately incorporated (p(x), v w (x), variable edge entropy, 
transition, roughness, etc.); (3) the higher order mean turbulence field models 
usually use a mean field model in the near-wall region (e.g., refs. 171 and 
172); (4) less computer time and storage are required, compared with mean turbu- 
lence field methods; (5) except for highly nonequilibrium flows, mean field meth- 
ods yield almost the same answers as mean turbulence field methods. (For exam- 
ple, see refs. 8, 24, and 173.) 


Details of the mean field closure can be conveniently discussed using 
sketch (a) where the boundary layer is shown subdivided into the usual three 
regions. The sublayer is the region nearest the wall. The no-slip and usual 
impervious wall boundary conditions ("wall discipline") require a wall damping 
expression as a modifier to whatever turbulence model is used in the other two 
regions. In the law-of-the-wall or fully turbulent region of the boundary 
layer, the turbulent motions are scaled as a function of y and experience indi- 
cates that the mixing length model is almost universally valid. The outer or 
wake region can be strongly influenced by "history" or relaxation effects and 
here the turbulent motions are scaled as a function of 6. 
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OUTER OR WAKE REGION 


LAW-OF-THE-WALL REGION 


SUBLAYER 

Sketch (a) 

Starting in the law-of-the-wall region, the most commonly used mean field 
closure model is the mixing length 

-u * v ' = l 2 3u 3u ( 37 ) 

By By 

For i a y (law-of-the-wall region) 

— u * v * = K 2 y 2 ( 38 ) 

By By 

where K is the Prandtl or Von Karman constant, which is approximately equal to 
0.4. 


In the outer, or wake region, there are two expressions commonly used, a 
mixing length expression 


-u' v' 


§ 2 A \ 2 3u 3u 


By By 


and an eddy viscosity expression (ref. 1 ) 


-u'v' = au e 6i — 


where, from computational experience (ref. 28 ), 6 ^ must be used instead of 6 * 

for the compressible case. 

For the sublayer, or wall damping region, several expressions are available 
(e.g., refs. 34, 40, 174, and 175), but results in reference 169 indicate that, 
at least for some cases, most of these expressions yield very similar results. 
The most commonly used wall damping expression was developed by Van Driest 
(ref. 174) and is an exponential damping upon the mixing length i(l - e~^^) 
where 


(A + = 26 for dp/dx and v w =» 0) 
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Therefore, in the wall region the mixing length model becomes 


-u' v' 


( Ky ) 2 ( 1 


e-y/A ) 2 


9u 9u 
3y 3y 


( 42 ) 


In addition to modeling the Reynolds stress term in equation (2), solutions 
of the compressible turbulent boundary-layer problem also require a model for 
the Reynolds heating term (eq. (4)). This is generally handled through use of a 
turbulent conductivity expression 

-vMt 7 = K t — (43) 

3y 


which then allows the definition of a turbulent Prandtl number 


K t = £ (44) 

»Pr,t 

Once Np r t is known, the model for -u'v' can be used in equation (44) to 
determine’ K t . 

Therefore, in the usual mean field closure models, there are four closure 
"constants" which must be determined or specified before a solution can be 
obtained: A + for the wall damping region, K for the law-of-the-wall region, 

( i / 6 ) mqy or a for the wake or outer region, and Np r?t . These constants will 
now be examined in some detail, with particular emphasis upon their possible 
variation due to (1) compressibility, (2) low Reynolds number, (3) wall blowing, 
(4) pressure gradient, and (5) roughness. 


Wall Damping Constant 

This is the most variable of the four "constants." Available data, 
although not definitive, show no discernible effect of Mach number upon A + 

(e.g., refs. 176 and 177). However, when A is computed from A + for compres- 
sible flow (A = A + v/\/t/p) , V, T, and p can be evaluated at several places 
(at the wall, locally as a function of y, and at the edge of the sublayer). 
Computational experience indicates evaluation at the wall is slightly better 
(refs. 53, 176, 177, and 178). 

There is no discernible influence of low Reynolds number upon A + , or at 
least none has yet been identified (e.g., ref. 179). The influence of wall blow- 
ing upon A + , however, is quite large, and there exist several ways of account- 
ing for this effect. One of the first procedures was a simple correlation plot 
of A + and the wall injection similarity parameter 2F/Cp (ref. 53). The data 
for this plot (fig. 26) were obtained from an examination of low-speed blowing 
data. By analogy to the F = 0 case, the A + -value was taken as the y + -location 
where the data were faired into a typical law-of-the-wall variation (i.e., the 
outer edge of the sublayer-buffer region). This variation (fig. 26) worked 
quite well (ref. 53), and most of the other models tend to give similar results 
(ref. 180). The changes in A + with wall mass transfer are quite large and 
some modification to A + must be included in order to obtain a reasonable C f 
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prediction. Probably the most general method of adjusting A + for blowing is 
given by Launder (ref. 171), who suggests A + = 26/(t/t w )'*7 (i.e., since 

x/x w = f (y) , A + = f (y) ) . 


The variation of A + with pressure gradient is also quite large. Cebeci 
and Smith (ref. 12) give a correction procedure, as does Kays (ref. 1 8 1 ) . (See 
also ref. 19.) Again, probably the simplest approach is to use the Launder 
expression A + = 26 /(t/t w ) 1-7 which accounts quite well for the effects of both 
mass injection and pressure gradient. The basic problem with the variation of 
A + with blowing and pressure gradient is indicated in sketch (b). As shown in 



the sketch, the wall damping occurs in a region where x 4 Constant, and the 
various approaches involve either letting A + n f(8x/3y) (ref. 179) or using 
the local x(y) to correct the A + or x w value; that is (from ref. 171, for 
example) , 


— 26 

(T/T w ) n 


( 45 ) 


Another parameter which has a large effect upon A + is roughness. 
McDonald and Fish (ref. 182) give a correction term for roughness as follows: 


Wall damping with roughness a (Wall damping without roughness) 


+ 


k + 

30(y + + 1) 



(46) 


where 



v 
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u T is based on t(y) and k is the roughness height. It should be noted that 
this damping term (with roughness) could be greater than 1. 

The influence of moderate acceleration (approaching laminar ization, involv- 
ing significant alteration of the sublayer) upon A + was investigated by 
Launder and Jones (ref. 183). They proposed the following expression: 

A + = 11 + 7900L (L > 0.0019) (47) 

where the pressure gradient parameter L is 

L = c; 2/3 _V_ 

u e 2 dx 

This expression has not yet received detailed consideration for compressible 
flows . 


Prandtl Wall Constant 

From available data, there is no appreciable effect upon K of either com- 
pressibility (general computational experience, ref. 184), low Reynolds number 
(refs. 179, 185, and 186), wall blowing (refs. 184 and 185), or surface rough- 
ness (refs. 45 (pp. 396-398) and 187). There is, however, a moderate influence 
of pressure gradient upon K, as shown in figure 27 (taken from ref. 188). From 
the work of Lewis, Gran, and Kubota (ref. 139), the Clauser pressure gradient 
parameter 8 r (and therefore fig. 27) is applicable to compressible flow if 6* 
is taken as 6*. Except for apparent pressure gradient effects (results similar 
to those in fig. 27 were also observed in ref. 30 (pp. 10-1 - 10-13) for super- 
sonic flow) , the wall slope K is the most constant of the mean field closure 
’’constants ." 


Outer or Wake Constant 

The present section discusses the behavior of (\/5) max (eq. (39)) instead 
of a; the behavior of a (eq. (40)) is quite similar, and in most instances 
the two expressions (eqs. (39) and (40)) give similar results. Figure 28 pro- 
vides a typical comparison between the use of equations (39) and (40) for calcu- 
lation of surface shear. For this low hypersonic case with moderate wall cool- 
ing, the mixing length is in slightly better agreement with the data. From the 
classic work of Maise and McDonald (ref. 79), there is no appreciable Mach num- 
ber effect upon U/<5) m ax* There is, however, a large low Reynolds number 
effect. Recent research (ref. 189) has strongly indicated that this low- 
Reynolds-number effect is caused by the residue of the boundary-layer transition 
process . 

One method of correlating the low-Reynolds-number effect is with a scaling 
parameter 6 + , where 

5+ _ gu T,max (48) 
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In references 176, 186, and 1 89 this parameter successfully correlated the 
low Reynolds number effect over a wide Mach number range. For the low-speed 
case a parameter Rg is sometimes used (ref. 175) when dp/dx = 0. Working 
with this expression for <5 + , one can show that 6 + a e/y or a turbulence 
Reynolds number as follows: 

*+ . 6 Vt # (4 9 ) 

W 

From the Townsend or Bradshaw assumption of the relationship between e and T, 

p u'v' « e p « T (50) 

Therefore, 

6+ « 6 P ^ (51) 

y 

From the Prandtl model for e 

£ “ l P \fe (52) 

and identifying i with 6 one obtains 

5+ cc £ (53) 

y 

The extent of the low-Reynolds-number amplification is indicated in fig- 
ure 29(a) (from ref. 189) where ( l /^)max was derived from experimental veloc- 
ity profiles for flat plates, cones, and cylinders. Increases of a factor of 2 
or more in (i/6) m ax above the usu al levels (0.07 to 0.09) for values of 5 + 
near 100 are observed. Since u'v 1 ~ (i/^max* this represents an increase of a 
factor of 4 or more in turbulent shear. The evidence that this increase is a 
function of distance downstream of transition is given by comparison of fig- 
ure 29(a) with figure 29(b) (also from ref. 189) where (l-/6) max actually 
decreases with decreasing 6 + for the nozzle wall case (where transition gener- 
ally occurs far upstream in the settling chamber). The correlation for inclu- 
sion of low Reynolds number effects should therefore be a function of both 6 + 
(or e/y) and Ax/6, the number of boundary-layer thicknesses downstream of the 
end of transition. (According to ref. 190, a value of Ax/6 of approximately 
30 to 50 is needed to "wash out" the low Reynolds number effect.) 

As Mach number increases, edge Reynolds numbers (such as R e x > R e 0 > 
etc.) become proportionally much larger than Reynolds numbers based on wall con- 
ditions (such as 6 + , etc.). This is due to the large difference, at least for 
wind-tunnel conditions, between the wall and boundary-layer-edge temperatures. 
Therefore, as shown quantitatively in figure 30, a given value of 6 + (which 
seems to correlate the transition-induced low Reynolds number effects) can cor- 
respond to a very large value of R e x . This implies that the low Reynolds num- 
ber amplification can effect a rather large portion of the boundary-layer flow 
on a hypersonic vehicle. 
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The outer or wake constant ( \ / 6 ) m a x is a] - so influenced by wall blowing 
and pressure gradient. From the data shown in figure 31 (from ref. 53) the 
level of (i/5) ma x decreases as the profile becomes less full and d6/dx 
increases. Reference 1 86 contains an excellent explanation for this effect, 
which involves the problem of flow history effects. A detailed discussion of 
this problem is postponed until later when nonequilibrium flows are considered. 
As a final point, for the outer region model there seem to be only moderate 
increases of (^/5)max with roughness (ref. 191). 


Turbulent Prandtl Number 

When deal ing with Np r £, there is a basi c qu est ion of wh ich term to model, 
v'H' or v'h' . Since H' = h 1 + u u* ,_then v'H' ~ v'h' + u u'v' and there- 
fore Np r _-r = e/iop (where v'H' = <p(3H/3y)) cannot equal Np r>t = e/< t (where 

v'h' = K^tSh/Sy)), except for M « 0. This question is discussed in some detail 
in references 53 and 192. The results of reference 177 clearly show that Np r> t 
is much more invariant with Mach number than Np r p. Typical Np r t distribu- 
tions inferred from profile data are shown in figure 32 for low-speed and super- 

sonic flow. Figure 32 shows data from references 181 and 193 to 196. An 
expanded scale is used, and the scatter, except near the wall, is typically 
+10 percent about a hypothetical mean line through the data. The scatter near 
the wall is probably caused by inaccurate mean profile data in that region. The 
data from reference 177 showing the problems with Np r p and Np r> ^ are shown 
in figure 33. The Np P ^ data at M = 7.2 agree fairly well with the low- 
speed data in figure 32, indicating that there may not be a strong Mach number 
effect on Np r> f 

The evidence for a possible influence of low Reynolds number upon Np r> t 
is primarily circumstantial, but indicates little or no influence (at least for 
air) . Figure 3^ summarizes some of the more accurate Reynolds analogy data 
(from refs. 197 to 201) for the low Reynolds number case. Since Np r ^ is a 
measure of the relationship between turbulent shear and turbulent heating, any 
large changes in Np r t might be expected to affect the Reynolds analogy fac- 
tor. The available data do not seem to exhibit a strong trend with 6 + , at 
least over the limited range for which data are available. 

A second piece of circumstantial evidence as to the effect of low Reynolds 
number on Np r j. is indicated in figure 35 (taken from ref. 12). The curves 
shown in this figure are the result of purely theoretical calculations based 
upon an assumed model of turbulence, but they do indicate only a very small 
effect of e/y upon Np r t for the Np r 0.7 case. However, if the low 
Reynolds number effect comes from the persistence of transitional flow struc- 
tures, this evidence is not really generic. In the absence of further data no 
firm conclusion can be drawn concerning the influence of low Reynolds number on 
Np r>t except that no large effect has yet surfaced. 

The variations of Np r t with y + (fig. 32) could also be interpreted as 
indicating a variation of ftp r t with e/y, but this influence of low Reynolds 
number near the wall results in reduced shear (e.g., ref. 202) as opposed to the 
influence of low Reynolds number in transitional flow structures, which results 
in increased shear. Therefore, as stated previously, e/y (or 6 + or y + ) is 
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not sufficient to adequately identify or correlate the various low Reynolds num- 
ber influences. These low Reynolds number problems obviously require further 
definitive research, especially in regard to Np F)t . 

The data from references 196, 203, and 204 indicate only a small influence 
of wall blowing on Np r In reference 204, data with near equilibrium pres- 
sure gradients indicate that Np r ^ values do become somewhat lower in the law- 
of-the-wall region for the positive dp/dx case. 


Intermittency Normal to the Wall 

Another ingredient often used in mean field closure procedures is a normal 
intermittency factor Ty, which is generally defined as the percentage of time 
that the flow is turbulent. This factor accounts for the superlayer or intermit- 
tent nature of the outer region in turbulent boundary layers. For high Reynolds 
number calculations ((e/y) max > 100 to 200) the inclusion of T y is generally 
a second-order effect (refs. 45 (pp. 366-374), 53, and 205); however, Ty can 
become important for (e/y) max < 100 (ref. 53). 

Typical distributions of Ty are shown in figure 36 (from ref. 25) for 
the low-speed and hypersonic cases. These data indicate that a fully turbulent 
flow occurs much farther out in the boundary layer for the high Mach number 
case, and that, even for the. simple dp/dx ~ 0 case, T y is a function of M. 

In addition to this Mach number effect, the data of Fiedler and Head (ref. 206) 
show that Ty can be a strong function of pressure gradient. The influence 
of roughness on Ty is evidently fairly small (ref. 207). 


APPLICATIONS OF MEAN FIELD CLOSURE METHODS 

Particularly during the last 8 years the literature has been rife with com- 
parison between theory and data for mean field prediction methods. These compar- 
isons indicate that mean field closure approaches yield predictions which are 
quite accurate over a wide range of conditions. (See especially refs. 12, 28, 
and 30.) In this section, comparisons of prediction and data for the more usual 
cases (such as dp/dx = 0, wall blowing, or equilibrium flows) will not be con- 
sidered, but instead, the more unusual situations, such as a diagnosis of nozzle 
wall boundary layers and transitional flows, among others, will be treated. 


The Case of the Nozzle Wall Turbulent Boundary Layer 

Thick, fully developed turbulent boundary layers have been difficult to 
obtain on models at high Mach number as a result of ( 1 ) the increase in transi- 
tion Reynolds number with Mach number (e.g., ref. 208), (2) the notorious diffi- 
culty in tripping high Mach number boundary layers, and (3) the general decrease 
in test section size of hypersonic facilities compared with the low-speed facil- 
ities. Therefore, a great deal of the hypersonic turbulent boundary-layer pro- 
file data were taken in nozzle wall boundary layers, which are readily available 
and usually of the order of 5 to 50 cm thick. 
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However, as can be readily seen with the help of sketch (c), there are sev- 
eral problems with these nozzle wall boundary- layer data. First of all, the 
measuring stations usually cover a small range of Ax/ 6 (data taken in test sec- 
tion region only with 6 quite large). These data are therefore not very satis- 
factory as test cases for the mean field methods, as these procedures are of the 
parabolic type and the solution proceeds in the streamwise direction. The data 
are usually taken at a specified x-location for a range of unit Reynolds number, 
rather than over a respectable range of Ax/ 6 at a given unit Reynolds number. 



USUAL MEASUREMENT 
REGION 

Sketch (c) 


The second problem involves so-called '’history" effects. Although the 
flow may locally inhabit a region where dT w /dx = 0 and dp/dx □ 0, just a few 
boundary- layer thicknesses upstream the gradients in these quantities are often 
quite large and there is a question of how fast the wall boundary layer loses 
the "memory" of these gradients. Of course, it is possible to make measurements 
at several stations along the nozzle, starting with a measured profile near the 
high gradient (large dp/dx) region (such as in ref. 47). This produces a quite 
interesting but theoretically undemanding test case (boundary-layer recovery 
from a favorable pressure gradient). 

An additional problem in the nozzle wall data is that of low Reynolds num- 
ber similitude. This problem concerns the early transition (usually in the set- 
tling chamber) and subsequent early "loss of memory" of the transitional flow 
structures. This problem was covered previously in connection with figure 29 
and indicates that the nozzle wall data for low values of 6 + may be applicable 
only to a limited class of applied flows, such as where roughness induces early 
transition on the nose of a vehicle before the boundary layer undergoes the 
expansion to afterbody flow conditions. 

Several investigators have studied these various nozzle wall problems (such 
as refs. 176 and 209 to 212). Figure 37 indicates typical total temperature and 
velocity relationships obtained on flat plates, cones, and cylinders (small 
dp/dx and T w (x) history), and figure 38 illustrates the nozzle wall case. 
These figures (taken from ref. 208) are now quite old (made up approximately 
8 years ago); however, the differences between the flat plate and nozzle wall 
data (that is, the differences between a linear variation of T^ and u as 
opposed to a quadratic variation) have also been consistently observed in the 
more recent data. The development of the quadratic variation is clearly seen in 
figure 39 (taken from ref. 212). 
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In summary, the current status of the nozzle wall data indicates that for 
M > 5, an adiabatic wall, and large values of 6 + (greater than 2000), the noz- 
zle wall profile data in the test section region are fairly typical (in approxi- 
mate agreement with the flat plate case). (See ref. 190.) For large values of 
S + , M > 5, and a nonadiabatic wall, the velocity profiles are still fairly 
typical (ref. 190), but the profiles are evidently out of equilibrium 

(ref. 209). For low values of 6 + (without relaminarization downstream of the 
throat), the profiles are evidently correct only for practical cases where tran- 
sition is far upstream. 


Influence of Wall Blowing Upon Skin Friction 

The conventional method of representing the influence of porous wall injec- 
tion upon skin friction is with a similarity plot of Cf/Cf 0 against 2F/Cf 0 . 
An exhaustive review article by Jeromin (ref. 33) presents the available data in 
these coordinates. (See fig. 1 of ref. 33-) Jeromin's plot indicates a strong 
Mach number effect upon skin friction reduction due to blowing. Several predic- 
tors (refs. 213 to 215) have examined the question of a Mach number effect in 
these coordinates, two of these since the publication of Jeromin’s review. 

Reeves (ref. 30 (pp. 6-1 - 6-A2-2)) calculated a weak Mach number effect, but his 
results indicated separate influences of Reynolds number and the ratio of wall 
temperature to total temperature, which could account for some of the apparent 
Mach number effects. Squire and Verma's (ref. 213) calculations (fig. 28 of 
ref. 213), in which a conventional mean field turbulence modeling procedure was 
employed, indicated very little influence of either Reynolds number, Mach num- 
ber, or T w /T£, at least for R e q > 8000 (the lower limit for their calcula- 
tions). This result is similar to conclusions reached in the early work of 

Rubesin (ref. 215). Landis (ref. 214) attributed most of Jeromin's Mach number 

effect to the influence of T w /T^. 

The data which indicate the strongest Mach number effect are those of 

Danberg (ref. 216). When the low Reynolds number effects are included in the 

analysis of Danberg's data (value of Cf 0 redefined), the final result for 
the available Mach number range (data ancl theory) is shown in figure 40 (from 
ref. 217). There is apparently very little Mach number effect, and the data are 
correctly predicted by the mean field closure methods. There is insufficient 
experimental data to determine any possible effect of the ratio of wall tempera- 
ture to total temperature. 


Transitional Flow Calculations 


The now classical approach to calculation of the transit iona l flow region 
between laminar and turbulent flow is to multiply the usual u'v' model by a 
streamwise intermittency factor T x to account for the increasing presence of 
turbulent bursts (e.g., refs. 142, 153, and 164). Using the definitions of x^ r 
and *tr,end indicated in sketch (d), the expression usually used for the inter- 
mittency ’ factor is (from ref. 2 1 8 ) 
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where 


(54) 


Ax _ x tr,end ~ x tr 

2.96 



Sketch (d) 


To apply this expression to a problem or a set of data, the x tr and 
x tr end v al ues must be known. Typical values of the ratio x^ r end /x tr (° r 
the’ ratio R^ r end^tr’ which is the same thing for dp/dx = 0) are shown in fig- 
ure 41 (from ref. 219). A nominal value of 2 for the ratio of x^ r en( j and 
xt r is seen to be a reasonable approximation over a large Mach number range. 

(See also ref. 201.) Therefore, the problem is reduced to a specification of 
xt r , a discussion of which is beyond the scope of the present report. (See, for 
example, ref. 220.) A comparison between equation (54) (which gives low-speed 
values of T x ) and hypersonic measurements of T x is shown in figure 42 
(from ref. 59). The high-speed data are in relatively good agreement with the 
curve from equation (54), indicating that the latter may be used with some 
confidence. 

For accurate transitional flow calculation in compressible flow, the 
precurser effects (discussed previously) should be included. An example of 
the increased agreement obtained when these precurser and low Reynolds number 
effects are included is given in figure 43 (from ref. 58). 


Additional "Nontypical" Applications of Mean Field Closure Methods 

Variable edge entropy .- This problem is caused by the mass "swallowed" or 
entrained by the boundary layer along the afterbody of slightly blunted vehi- 
cles. (See sketch (e).) As the boundary layer swallows the high entropy stream- 
lines, the edge properties possess a variable entropy condition. The actual 
edge conditions are determined by equating mass flow in the boundary layer at a 
given body station m^ to the mass flow in an entering stream tube ^shock • 

This problem was treated in reference 221, using a mean field closure aproach. 
Figure 44 (from ref. 221) indicates the better agreement resulting from consider- 
ation of variable entropy. 

Transverse curvature influence .- Cebeci and Smith (ref. 12) have treated 
this case quite well, using a mean field approach. The problem arises when 
6/r c ~ 0(1), where r c is the local body transverse radius of curvature. 
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According to reference 12, the outer region formulation is unchanged. The inner 
and wall damping regions are altered to 



Using this expression Cebeci and Smith obtained good agreement with data on a 
0.06l-cm-diaraeter needle at M n 5.8. 

OO 

Longitudinal curvature influence .- The definitive review work in this area 
is a fairly recent report by Bradshaw (ref. 222). The usual effect is for con- 
vex curvature to reduce the Reynolds stress and entrainment and for concave cur- 
vature to increase the shear over the no-curvature level for the same pressure 
gradient. There are evidently four distinct effects of longitudinal curvature 
upon compressible turbulent boundary- layer physics and calculations. The first 
of these effects is the influence of additional longitudinal curvature terms in 
the mean flow equations. These extra terms were included in the calculation 
methods of references 167 and 223 and calculations in reference 224 using the 
code of reference 1 67 ; the results of these calculations indicate better agree- 
ment for a flow with mild concave curvature when these terms are included. 

The second effect is unique to compressible flows, in that longitudinal 
curvature can induce a large pressure gradient in the boundary layer. This 
influences the flow even beyond the boundary layer and creates a nonuniform free 
stream (i.e., induces u e (y) for y > 6). This occurs primarily as a result of 
the Prandtl-Meyer type of flow turning which is a typical feature of compressi- 
ble flows. Therefore, the expression 3p/3y =0 is no longer true, and the cor- 
rect p(x,y) behavior must be included in boundary-layer prediction procedures; 
that is, one must have the correct value of Piocal to com P ute a density, and 
3p/3x is now a function of y. The codes of references 1 67 and 223 both 
include this capability. 
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The third effect of longitudinal curvature is the influence on the turbu- 
lent structure itself. Bradshaw (see ref. 222) has given a correction factor of 
the form 


i = (l\ A - B ( 56 ) 

6 U) t \ du/dyj 

where 

B « 0(10) 

This term is generally applied in an average sense to the outer region flow 
using the wall curvature for r c and provides a first-order correction for mod- 
erate curvature (herein defined as convex radius of curvature large enough so 
that laminarization does not occur and concave radius of curvature large enough 
so that Goertler vortices are not important and separation does not occur). 

The fourth effect of longitudinal curvture is illustrated by experimental 
data, which indicate that for 6/r c > 0.005 (which is not very large), concave 
curvature can generate steady, longitudinal Goertler vortices embedded in the 
outer portion of the turbulent boundary layer (ref. 225). The presence of 
Goertler vortices converts a readily solvable two-dimensional boundary- layer 
problem into a more complex three-dimensional turbulence problem of the 
parabolic-elliptic type (recirculation in the crossplane). Much more research 
is needed before quantitative predictions about the latter type of flow can be 
made. 


Calculations for adverse pressure gradient in compressible flow .- There is 
currently some controversy concerning the capacity of turbulent boundary-layer 
calculation methods to compute adverse pressure gradient flows for the compressi- 
ble case. Bradshaw, in reference 226, partly on the basis of disagreements 
between his procedure and data, suggests that a mean dilitation correction is 
needed before satisfactory results can be obtained. Part of the problem in this 
area is due to the fact that the early (pre-1969) data for adverse pressure gra- 
dient compressible turbulent flows were taken on bodies where positive values of 
dp/dx were induced by longitudinal curvature, thereby introducing all the possi- 
ble problems just discussed. However, Bradshaw (in ref. 226) considers the 
newer data, where the waves causing the adverse pressure gradient are impressed 
upon a flat surface. 

One example of an adverse pressure gradient calculation (on a flat surface) 
is given in figure 45, and the external Mach number distribution for this flow 
is given in figure 46. The data are from reference 139 and the mean field calcu- 
lation method is taken from reference 12. In this case the comparison is quite 
good. A check calculation made by the present authors using the method of ref- 
erence 53 yielded similar agreement with the data. On the basis of these 
results, the question of whether Bradshaw’s dilitation correction factor 
(ref. 226) is really required to compute adverse pressure gradient compressible 
flow is still open. 


NONEQUILIBRIUM AND MEAN TURBULENT FIELD CLOSURE 


Physical Problem 


Primarily as a result of the large difference in scales between the inner 
(near wall) and outer portion of turbulent boundary layers, the inner portion 
(with scales of the order of y) reacts much faster to a change in boundary con- 
ditions than the outer region (with scales of the order of 6). This is shown 
schematically in sketch (f). The value of Ax/6 necessary for the outer region 
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Sketch (f) 


flow to "relax" and equilibrate with the new boundary conditions is of the order 
of 30 to 90. 

Therefore, as discussed earlier, for rapid chan ges in boundary conditions a 
rate equation is needed for the computation of u'v' in the outer region (and 
perhaps for A + as well; see ref. 1 8 1 ) . The turbulent field cannot instanta- 
neously follow or track rapid changes in the mean flow. A nonequilibrium situa- 
tion can be set up by the sudden removal or application of a pressure gradient, 
wall injection or suction, wall temperature gradient, or wall roughness. Larger 
changes, such as those caused by shock-wave impingement and large obstacles or 
steps, may cause the flow to violate the boundary- laye r assumptions. In addi- 
tion to the computation of nonequilibrium flows, u'v' rate equations and mean 
turbulence field models are also necessary for including free-stream disturbance 
boundary conditions and for calculating the various turbulent second-order 
correlations. 


There is a zeroth-order nonequilibrium modeling which involves a rate equa- 
tion for the mixing length or outer region constant. This equation is usually 
of the form 

^ Hx = ^ ee Q u ilibrium ” (57) 

Simple expressions or rate equations of this general form were used in refer- 
ences 30 (pp. 6-1 - 6-A2-2 and 29-1 - 29-10), 45 (pp. 375-383), 227, and 228 and 
do tend to give at least the correct qualitative behavior, the problem being 
that the decay constant can change with the particular flow involved; that is, 
the equation has the correct form but the absolute rate of return to equilibrium 
is uncertain. Therefore, the balance of this section is devoted to methods 
based upon rate equations derived from the Navier-Stokes equations (second-order 
correlations or Reynolds stress transport equations). 

Second-order equation closure methods (or mean turbulence field methods) 
have been applied rather extensively to low-speed turbulent boundary layers. 
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(See refs. 21, 172, and 229 for reviews of this work.) However, their applica- 
tion to compressible flows has been rather limited, perhaps as a result of the 
following: (1) relative success of the mean field methods, (2) dearth of accu- 

rate data for nonequilibrium compressible turbulent boundary-layer flows (i.e., 
very few interesting test cases), (3) questions concerning possible importance 
of p' and p' terms in these equations, and (4) lack of sufficient data on 
the fluctuation field for the compressible case (for modeling purposes). 


Equations and Modeling for Compressible Nonequilibrium Flows 

The Reynolds stress transport equations (or Rj^j equations) are generally 
obtained by (1) multiplying the Navier-Stokes equations for by uj, (2) mul- 

tiplying the Navier-Stokes equations for uj by u it (3) adding these two equa- 
tions, (4) substituting for instantaneous values the usual mean plus fluctuation 
portions and Reynolds averaging, (5) substracting out the kinetic energy equa- 
tion for the mean flow, and (6) imposing the usual boundary-l ayer order of mag- 
nitude analysis. The result of all this is an eq uation for u|uj. From this 
equation o ne c an obtain (1) an equation for uj 2 (i = j = i } 2, 3), (2) an equa- 
tion for u'v* (i = 1, j = 2), and (3) an equation for e = u' 2 + v' 2 + w’ 2 

(sum of equations for uj 2 ^. The general form of these second-order equations 
is usually the following: 

Convection = Production + Diffusion - Dissipation 

Details of the derivation of these equations are available, for example, in ref- 
erences 36, 42, 148, 159, and 230 for the compressible turbulent kinetic energy 
equation and in references 43 and 156 for more complete sets of the second-order 
compressible equations. Reference 12 represents probably one of the best refer- 
ences for general discussion and derivation of the second-order equations for 
both low and high speeds. 

From computational experience, the critical portion of the development of 
a good nonequilibrium turbulent boundary- layer calculation procedure is an accu- 
rate equation for the length scale of turbulence (refs. 36 and 45, pp. 275-299). 
The results of the calculations for nonequilibrium flows are usually not 
extremely sensitive to the details of the modeling used in the second-order cor- 
relation equation, but the length scale determination must be nearly correct for 
a good prediction. 

There are three approaches usually used to derive a length scale equation: 

( 1 ) form an equation for the turbulence dissipation D and relate this to L 
through the usual model for D (D * e3/ 2 /L) (ref. 8); (2)' use an equation for 
the two-point correlations, which introduces a length scale quite naturally; and 
(3) use the turbulent vorticity equation and relate to L, for example, through 
the expression L = (refs. 231 and 232). An alternate approach to the 

length scale problem is to use an algebraic relationship (such as a function of 
y/6, for example), which is assumed known (ref. 148). This latter method does 
not yield satisfactory results for truly nonequilibrium flows (if the length 
scale being solved for is the one used in modeling the Reynolds stress) . 
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Single-Equation Models 


Bradshaw approach .- This method (ref. 36) gives a second-order outer region 
solution based upon the conversion of the turbulent kinetic energy equation into 
an equation for u' v' by assuming that a-) = x/2pe = Constant. The only length 
scale which appears is in the dissipation term; this is a very important aspect 
of the method, as Huffman (ref. 233) has shown that this length scale L (for 
dissipation) is much more invariant than the usual mixing length (which is used 
to model turbulence production), as can be seen by comparing figures 47 and 48 
(from ref. 233). Bradshaw's a-| is also relatively invariant, even for nonequi- 
librium flows, as indicated in figures 49 (from ref. 233) and 50 (from ref. 234). 
Bradshaw, by using a-| , evidently sidesteps the problem of modeling the produc- 
tion term with a length scale and then having to derive and solve an equation 
for that scale. The length scale which Bradshaw uses (in the dissipati on t erm) 
is evidently fairly invariant, and therefore, his one-equation method (u'v' equa- 
tion plus algebraic length scale) should yield predictions which are as accurate 
as the more common two-equation models (e and i equations) . 

Method of Shamroth and McDonald and method of Chan .- These procedures 
(refs. 159 and 155) use the integral form of the turbulent kinetic energy equa- 
tion and Bradshaw modeling (a-| and L). The procedures are, at least approxi- 
mately, integral forms of the Bradshaw closure approach and seem to be fairly 
successful in predicting nonequilibrium flows. The method is inherently fast 
and relatively accurate and is a good choice for a simple nonequilibrium 
boundary-layer calculation procedure. 


Two-Equation (and More) Models for Compressible Flows 

Method of Wilcox and Alber .- This procedure (ref. 231) solves the turbulent 
kinetic energy equation and the fluctuating vorticit y eq uation. Mass-weighted 
averages are used along with the Prandtl model for u'v' ~ e^^i(3u/9y) . 

Method of Spalding and Gibson .- This method (ref. 232) is quite similar to 
that of reference 231 in basic concept (e, w equations used) but quite differ- 
ent in detail and application. Neither of these procedures has yet seen much 
application in highly nonequilibrium compressible flows. 

Method of Donaldson and Sullivan .- This method (re f . 156 ) s olv es the c om- 
plete set of second-order boundary-layer equations (u' 2 , v' 2 , w' 2 , u'v', 

T'T' , u'T' , and v'T' ) and therefore does not have to use a turbulent Prandtl 
number. However, most of the p' type terms were neglected. The latest ver- 
sion of this method (in publication) includes a length scale equation. The type 
of detailed information obtainable from this type of closure is shown in fig- 
ures 51 and 52. 

At the present time there are insufficient detailed data available to 
develop "calibrated" nonequilibrium (mean turbulence field) closure methods for 
compressible turbulent boundary layers. 
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CALCULATION OF THREE-DIMENSIONAL COMPRESSIBLE 


TURBULENT BOUNDARY LAYERS 
Status 

A general method for solving three-dimensional compressible turbulent 
boundary-layer flows is not currently available; however, a limited number of 
problems have been solved using special techniques for each case. Blottner 
(ref. 235) presents a review of specialized and approximate solutions of the 
laminar three-dimensional boundary- layer equations. General three-dimensional 
boundary-layer computer codes require the following: ( 1 ) an efficient and accu- 

rate numerical procedure, including the required logic to automatically change 
the difference molecule as a function of local flow conditions, (2) a general- 
ized curvilinear coordinate system convenient for the design engineer, (3) accu- 
rate turbulence models for the Reynolds stresses, and (4) a three-dimensional 
inviscid flow field solution which is compatible with the boundary-layer coordi- 
nate system. 

A number of general numerical procedures and computer codes have been 
developed for two-dimensional and axisymmetric turbulent boundary- layer flows 
(e.g., see refs. 53, 142, 1 64 , 202, and 236 to 242). This success was extended 
to a particular class of three-dimensional flows designated as quasi-two- 
dimensional (retains computational advantage of two independent variables while 
allowing large cross-flow components) by Hunt, Bushnell, and Beckwith (ref. 243). 
Numerical solutions for general three-dimensional flows did not materalize until 
after 1972 because of computer systems limitations (storage and processing speed) 
and the lack of accurate three-dimensional inviscid flow field solutions required 
for boundary-layer-edge input conditions. 

Since 1972 an intensive research program on three-dimensional boundary- 
layer flows has been under way at a number of research centers. This effort is 
stimulated by increasing awareness of potential savings in cost and man-hours 
through numerical simulation of complex flows; often these flows cannot be simu- 
lated in ground test facilities (such as, for example, real-gas boundary-layer 
flow over the space shuttle). Substantial progress has been made as a result of 
the increasing availability of large-storage, high-speed computer systems and 
the maturing of accurate numerical procedures and computer codes for solving the 
three-dimensional inviscid equations for complex configurations. (See refs. 244 
and 245.) Some of the more important areas of three-dimensional boundary-layer 
research are as follows: coordinate systems and transformations, numerical solu- 

tion techniques, turbulence modeling, geometry, initial conditions, inviscid 
boundary conditions, inflow lines, regions of influence and dependence, and 
numerical optimization for perfect-gas flows as opposed to real-gas flows. 

These specific areas will be discussed in detail in subsequent sections of the 
present review. 

The primary purpose of this section is to discuss in detail problems associ- 
ated with development of three-dimensional compressible turbulent boundary-layer 
codes for application to complex aeronautical and aerospace vehicles and to indi- 
cate progress to date in achieving this goal. Comparisons of numerical results 
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with experimental data for several operational three-dimensional boundary-layer 
computer codes will also be presented. 

Problems associated with obtaining accurate three-dimensional flow field 
results will not be discussed. In reality, the prediction of the inviscid flow 
field with the accuracy required presents a challenge as difficult as solving 
the three-dimensional compressible turbulent boundary-layer equations; however, 
the inviscid flow field input will be assumed available in the present review. 
(See, for example, ref. 246 .) 


Boundary-Layer Equations 


An orthogonal curvilinear coordinate system is generally used to define the 
surface over which the boundary layer is flowing. The coordinate normal to the 
surface is X3 with X3 being zero on the surface. The lines x-| = Constant 
and X2 = Constant define the system of orthogonal coordinates on the surface. 
(See fig. 53 .) The square of the element of arc (ds) in the boundary layer is 

ds 2 = (h-| dx -]) 2 + (h2 dx2 ) 2 + (hg dX3 ) 2 (58) 

where h-| , h2, and I13 are metric coefficients (or scale factors): 

h-| = h-|(x-|,X 2 ) h 2 = h 2 (x-|,x 2 ) ( 59 ) 

The metric coefficient I13 is generally assumed to be unity (113 = 1) since the 
boundary layer is assumed to be thin. The governing equations (first order) can 
then be written as follows, where 113 is included for generality: 

Continuity 


_i_(h2h3pu-| ) + _i_(h-|h3pu2) 
3xi 3x2 


_^_(h-|h2PU3) = 0 


3 x 


( 60 ) 


x-) -momentum 
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X2-momentum 
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(62) 


Energy 


P U 1 3H + P^2 3H_ + P^3 9H_ = 1 <L 

h-j 3 xi h 2 3 x 2 hg 3 x^ I13 8x3 


3H 


Np r h3 8x3 


+ ]lh - _L\J 3 /u-! 2 + u 2 2 \ _ , H , 

\ N Prjh 3 3 x 3 \ 2 ) 5 


( 63 ) 


where 


PU3 = PU3 + p'u^ 

The geodesic curvatures of the surface coordinate lines are 

K, = - _L_ 

hih 2 3 x 2 


(64) 


k 2 = - _L_ 

h-|h 2 3x-j 


(65) 


The governing equations are completed with the perfect-gas equation of 
state (see refs. 245 and 247 for real -gas flows) 

p = pRT (66) 

a viscosity law 

U = f (T) (67) 

and suitable models or transport equations for the Reynolds stress terms appear- 
ing in the equations. 

Obvious differences between the equations for three-dimensional (eqs. (60) 
to (67)) and two-dimensional (eqs. (1) to (4)) boundary-layer flows are the 
appearance of the cross-flow terms in the continuity and streamwise momentum 
equations and the addition of a second momentum equation for the cross-flow 
direction. However, the most important difference between two- and three- 
dimensional equations is that the three-dimensional equations are hyperbolic 
rather than parabolic in coordinate planes parallel to the wall boundary. The 
hyperbolic character of the three-dimensional equations has been recognized by 
a number of authors (see refs. 248 to 257) and arises naturally from the varia- 
tion of the cross flow along the coordinate line normal to the wall boundary 
(x^ d Constant, x 2 = Constant). Streamlines originating at different points 
along this normal line diverge in the downstream direction. These streamlines 
originate at different upstream locations; consequently, wedge-shaped regions 
of influence and dependence are associated with three-dimensional viscous flow 
and extend upstreani and downstream from the computational point along the normal 
to the wall boundary. This characteristic of the three-dimensional boundary- 
layer equations is generally referred to in the literature as the influence 
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principle. (See ref. 248.) A schematic showing the regions of dependence and 
influence is presented in sketch (g). A more complete discussion is presented 
by Raetz (ref. 248), Kitchens, Gerber, and Sedney (refs. 254 and 257) and Wang 
(ref. 256 ). 



The required boundary conditions for the governing equations (eqs. (61) to 
( 63 )) are as follows: 

= 0 (wall) 


u-] = 0 

U2 = 0 


u 3 = u w (x 1 ,x 2 ) 

T = T w (xi,x 2 ) or (p\ = q w (xi ,x 2 ) > 

V x 3/w 

X 3 -*• 00 (edge) 

U 1 = u 1,e( x 1* x 2) 

u 2 = u 2,e( x 1> x 2) 

T = T e (x-| ,x 2 ) 


( 68 ) 


The governing equations for the inviscid flow at the outer edge (x^ -*■ °°) 
are evaluated from equations ( 61 ) to ( 63 ) as follows, where the static tempera- 
ture is utilized instead of the total enthalpy: 


x^] -momentum 

u 1 ,e 9 u 1 ,e + ^e 9 u 1 ,e 
h-| 9x 1 h 2 3x 2 


2 

u 1 ,e u 2 ,e K 2 + u 2 , e K 1 


1 3p 
Pe h 1 9x 1 


(69) 
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(70) 


X2-momentum 

U 1 ,e 9u 2,e + 9u 2,e 

hi 8x2 h£ 3x2 


2 

U 1 ,e u 2,e^1 + U 1 , 2^2 


_J 8p 

Pe h 2 9x 2 


Energy 

Pr n n f u 1,e 9T e + ^e\ = 8p_ + 8p_ (71) 

h-j 8x-| h2 8x2 / hi 3xi h2 8x2 

Variable entropy and vorticity can be treated by the system of equations; how- 
ever, its proper treatment is very complex for three-dimensional flows. The 
reader is referred to Mayne (ref. 255) and Popinski and Davis (ref. 258) for 
blunt cone flows and to Kendall et al. (ref. 245) for general configurations. 
The proper treatment of variable entropy is of primary importance for some 
applications such as the space shuttle. 


Two general approaches may be followed in evaluating the required edge con- 
ditions: (1) the pressure gradient terms required in equations (61) and (62) 

(dp/dx-L, i = 1, 2) can be obtained by direct solution of the Euler equation, or 
(2) the pressure p(xi,X2) can be specified together with the appropriate bound- 
ary and intial conditions and u^ e and T e obtained directly from the solu- 
tion of equations (69) to (71). the optimum approach is to develop simulta- 
neously an accurate three-dimensional inviscid computer code with a coordinate 
system compatible with that developed for the three-dimensional boundary-layer 
program; in any event, accurate and consistent edge conditions are required if 
accurate boundary- layer solutions are to be obtained. 

The sufficient conditions required to start the boundary-layer solution 
have been presented by Ting (ref. 259). In principle the system of equations 
can be solved numerically by marching parabolically away (downstream) from spec- 
ified initial data planes; however, the influence principle must be carefully 
treated. 


Coordinate Systems and Transformations 

Coordinate systems .- A number of papers have been published on the solution 
of three-dimensional boundary- layer flows; however, these papers have generally 
used specialized coordinate systems for each test problem. The coordinate sys- 
tem was generally chosen to simplify the governing equations for the particular 
flow of interest. The primary concern in the selection of a coordinate system 
is that it should allow the procedure to start from the initial conditions and 
proceed in a logical fashion over the entire surface without having to match 
together two or more distinct coordinate systems. 

The majority of the references available have used streamline coordinates. 
(See ref. 260, for example.) The streamline coordinate system is an orthogonal 
surface coordinate system formed by the projection of the inviscid streamlines 
and their orthogonal trajectories on the surface. The system has a number of 
distinct advantages; however, its calculation is a major numerical effort in 
itself and must unfortunately be repeated for each change in flow conditions 
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(angle-of-attack changes, for example). Another disadvantage is encountered in 
the design of wings and control surfaces where displacement surface effects must 
be accurately treated (for example, supercritical wing design). The coordinate 
system has to be recalculated for each iterative cycle between the three- 
dimensional inviscid and boundary-layer programs. This represents a major under- 
taking and would be prohibitively expensive in terms of computer processing time. 
Streamline coordinates also tend to diverge greatly in highly three-dimensional 
flow, which results in large truncation errors unless additional streamlines are 
introduced in order to retain a reasonable mesh point spacing in the direction 
normal to the streamline coordinate. Consequently, it should be realized that 
while streamline coordinates are optimum for many simple flows (for example, 
three-dimensional stagnation point flows), they do not appear at the present to 
be optimum for general three-dimensional flows. It is then advantageous to seri- 
ously consider the development of a body-oriented surface coordinate system that 
would completely eliminate the requirement of coordinate system recalculation 
for each change in flow condition. 

A body-oriented surface coordinate system could be made compatible with the 
inviscid three-dimensional computer program used to specify the required edge 
conditions. The only disadvantage of an orthogonal surface coordinate system is 
that the initial data lines cannot in general be made to coincide with the body 
coordinate lines even for simple geometry, such as a blunt nose cone at angle of 
attack. Blottner and Ellis (ref. 26 1) have presented an orthogonal surface coor- 
dinate system whose origin is at the stagnation point for analytic bodies at 
angle of attack. McGowan and Davis (ref. 262) have utilized an orthogonal sur- 
face coordinate system for sharp right circular and elliptical cones at angle 
of attack. However, these applications are for analytic bodies of revolution 
and not for general aerodynamic vehicles. The body-oriented coordinate system 
appears optimum from the viewpoints of geometry input and direct coupling with 
existing three-dimensional inviscid computer codes. The system completely 
avoids the time-consuming and difficult problem of streamline trajectory calcu- 
lation for each change in flow conditions; consequently, inviscid-viscous cou- 
pling procedures can be greatly simplified for problems where displacement sur- 
face effects are important. The difficulties associated with developing 
efficient procedures for numerically generating the coordinates and required 
metric coefficients for an orthogonal body-oriented surface coordinate system 
are more than compensated by program flexibility from the viewpoint of the engi- 
neer who must use the computer program as a design tool. 

The body-oriented surface coordinate system can be chosen to be either 
orthogonal or nonorthogonal . Regardless of whether the system is orthogonal or 
nonorthogonal , the X3~coordinate must be chosen normal to the wall boundary. 
Consequently, for the orthogonal system the selection of one coordinate x-| or 
X2 automatically fixes the remaining coordinate because of the orthogonality 
requirements. For example, if x-| is chosen to lie along the rays of a sharp 
right circular cone, then the X2~coordinate for x-| = Constant cuts the cone 
in circular elements; however, for sharp elliptical cones a warped cross-section 
element is formed. (See fig. 2 of ref. 262.) The nonorthogonal coordinate sys- 
tem appears desirable from a number of viewpoints: (1) It can be made to coin- 

cide with the location of the initial data planes; (2) it can be generated as 
easily as the orthogonal system; and (3) the difference grid can be developed to 
completely cover the computational region of interest - for example, swept, 
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tapered wings. (See fig. 2 of ref. 244.) Further, the nonorthogonal system 
results in only a minimal increase in the number of correction terms arising 
from the metric tensor (see ref. 245) provided the cross-flow diffusion terms 
are not considered (in other words, if only the classical three-dimensional 
boundary-layer equations are used). 

From the viewpoint of the design engineer it appears that a general compu- 
ter code for solving the three-dimensional boundary-layer equations should be 
developed in nonorthogonal, body-oriented surface coordinates. This system 
would also allow complete flexibility and could utilize streamline coordinates 
if so desired for special flow cases. This viewpoint is currently being incor- 
porated in the computer program under development by Cebeci et al. (ref. 244) 
for wing surfaces. 

Transformations . - The selection of suitable transformations for three- 
dimensional laminar or turbulent boundary- layer flow presents a problem since no 
one transformation can be considered general. The selection of primitive varia- 
bles introduces problems since the boundary layer is thin and grows at different 
rates in the x-|- and X 2 -coordinates; furthermore, the equations are singular 
at the tip of sharp bodies. When physical coordinates are used, the solution is 
extremely sensitive to the surface mesh point distribution; consequently, proce- 
dures developed in physical coordinates generally require either excessively 
small grid spacing distributions for the x-|- and X 2 ~coordinates or special 
treatment of the spatial derivative terms in the finite-difference equations if 
accurate results are to be obtained. This results directly in excessive compu- 
ter storage and/or processing time. Specific problems associated with physical 
coordinates are as follows: ( 1 ) excessive growth of the boundary layer in the 

computational domain requiring additional mesh points in the X3~coordinate as 
the solution proceeds downstream, ( 2 ) large gradients in the wall region for tur- 
bulent flows requiring either closely spaced mesh points in the normal direction 
or the inclusion of the law-of-the-wall relationship (not really satisfactory 
for three-dimensional flow), and ( 3 ) singularities such as at the tip or leading 
edge of sharp-edge bodies. The problem associated with the boundary-layer 
growth can be treated with suitable stretching and normalization as presented by 
Kendall et al. (ref. 245). 

Transformation to similarity variables has been shown to be useful for com- 
puting similar flows (see ref. 263 , for example); however, their usefulness is 
questioned by some authors (see, for example, ref. 245) for three-dimensional 
flows which are highly nonsimilar in character. The majority of the three- 
dimensional boundary-layer solution procedures now in existence use similarity 
type variables. (See ref. 245.) 

As in two-dimensional boundary layers the physics of turbulent flow pre- 
sents a numerical problem since wall gradients are large in comparison to laminar 
flow; the viscous sublayer generally requires a minimum of two to three mesh 
points. Generally a geometric-progression grid point spacing normal to the wall 
is employed; that is AX 3 ^ + -| = K AX 3 where K = 0(1.02). It has also been 
suggested that a logarithmic spacing &e utilized as indicated by the law of 
the wall. The main object, regardless of the transformation or stretching 
selected for the normal coordinate, is to minimize the total number of mesh 
points normal to the wall boundary in order to minimize the required computer 
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storage and processing time. It should also be noted that the x^-transformation 
generally complicates the Reynolds stress models and their inclusion in the com- 
puter code; this is one advantage of using primitive variables with suitable 
thickness normalization. (See ref. 245.) 

To date no general transformation has been developed for three-dimensional 
turbulent boundary-layer flows. Several transformations have been studied and 
are presented by Blottner (ref. 242). Some of these transformations will be dis- 
cussed in a subsequent section where three particular computer programs will be 
presented, with numerical results compared with experimental data. 


Turbulence Modeling 

Research on turbulent transport equation models and turbulence structure 
(e.g., refs. 30 (pp. D-1 - D-14 and B-1 - B-24) and 264) is very important since 
the goal of this research is to develop one set of equations with one set of con- 
stants or functional values which will adequately model all turbulent flows; how- 
ever, to date no transport equation procedure has been shown to be superior to 
simple eddy-viscosity models for two-dimensional, equilibrium, boundary-layer 
flows. (This is not the case for free shear, jet, and wake flows.) It is 
agreed by most investigators of turbulence modeling that eddy-viscosity models 
leave a great deal to be desired since the physics of the flow is in general 
neglected; however, from the des,ign engineers viewpoint, the eddy- viscosity 
concept works well for a wide class of two-dimensional turbulent boundary-layer 
flows. 

A number of rather difficult questions must be answered concerning the 
extension of the two-dimensional eddy-viscosity models to three-dimensional 
flow; for example, the cross-flow momentum equation requires the specification 
of an eddy- viscosity coefficient for the cross-flow direction. The simplest 
approximation is to assume that e X 2 = e x -| (isotropic eddy viscosity). This 
simple assumption implies that the turbulent stresses and mean flow rate of 
strain are always parallel . However , measurements of the complete turbulent 
stress tensor by Van den Berg et al. (ref. 265) indicate that the turbulent shear 
stress level in three-dimensional boundary-layer flow is lower than would be 
expected from extrapolation of two-dimensional theory; also, the local shear 
stress direction does not in general coincide with the direction of the velocity 
gradient. Experimental data indicate that the eddy viscosity should be lower 
in the cross-flow direction; that is, e x2 //e x1 < ^ • One might assume that 
e x2 = ae x1 where a < 1; however, this approximation is unsatisfactory (see 
ref. 266) unless some empirical functional relationship between e x i and e X 2 
can be determined which allows accurate numerical prediction of a wide class of 
three-dimensional boundary-layer flows. This approach may become feasible as 
accurate three-dimensional experimental data are obtained. Two-dimensional eddy- 
viscosity models were developed to their current level primarily as a result of 
the availability of a wealth of experimental data; the same data base must be 
provided for three-dimensional flows for successful models to be developed. 

In the present paper, closure through the solution of transport equations 
will not be discussed in further detail. The reader interested in this approach 
is referred to references 172 and 267 to 269 for examples. 
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Closure of equations (61) to (63) requires that an effective viscosity 
y e ff and thermal conductivity be expressed in terms of the mean flow 

variables. The Reynolds stress terms are defined as follows: 


“l — _ 8ui^\ 

-P U 1 U 3 = e x1 

8x3 

» r 3up 

-p U 2 U 3 = ex2 

-P ^H 7 = £0 |S- 

3x 3 J 


(72) 


The further assumption is made that e x i = e x2 = e and a turbulent Prandtl num- 
ber is defined which relates the effective conductivity and effective viscosity. 
The following relations can then be formulated: 


Mef f = y + jj rj 


(73) 


and 


K eff 


°P y ft + e N Pr A 
Np r \ ii N Pr . ft J 


m) 


The streamwise intermittency function T x (ref. 218) models the transitional 
region of flow and is a function of x-| and x 2 (0 ^ T x ^ 1; T x = 0 for lami- 
nar flow; T x = 1 for turbulent flow). The effect of pressure gradient on inter- 
mittency is presented in reference 270. For the present discussion it is 
assumed that the initiation and completion of the transitional flow region are 
specified; however, correlation relations could be directly incorporated into a 
computer code. 


The simplest eddy- viscosity approach (see fig. 54) is to assume that 
the eddy viscosity is a scalar function independent of coordinate direction 
(refs. 168 and 243). The following models are considered since they are used 
in most current three-dimensional, compressible, turbulent, boundary-layer 
programs : 

Single-layer model 



where 

11 = k 2 tanh/!H x 3 \ d (76) 

x 3,e \ k 2 x 3,e/ 
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1 - expf- ^3 


a+ m fly" 

\P/w\P /w 


1 - erf [k 5 x 3 - k6 


= (t*1 


2 \1/2 
?x2 


Two-layer model 


Inner law 


= P*2 


Outer law 


,2 ± + 


h 3 2 |_\ 3x 3 


where 


e outer = k 4P v e^*^y 


12 = icx^D 


(0 g x 3 $ x 3;C ) 


(x 3 , c < x 3 g x 3>e ) (82) 


/ 2 2 
V e = \ U 1 , e + u 2 , < 


r x 3,. 

" -'0 


( Ul 2 + u 2 2 ) 


The point where the inner and outer laws are matched x 3>c is obtained from the 
continuity of eddy viscosity; that is, ^inner = e outer- ( See sketch (h).) No 
attempt is generally made to assure continuous derivatives for the two laws; how- 
ever, the single law is continuous and appears to be as satisfactory as the more 
complex two-layer relationship (ref. 271). Three-layer models have been pro- 
posed by Pletcher (ref. 158), but there does not appear to be any advantage to 
the multilayer models above the single-layer mixing length model (eq. (75)) or 
the two-layer model (eqs. (81) and (82)). The normal intermittency factor 
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Sketch (h) 


(eq. (79)) proposed by Klebanoff (ref. 272) is generally approximated by the 
following relationship: 


-1 

(86) 


For high Reynolds number flows the constants appearing in equations (75) to 

(83) are generally assigned the following values: k^ = 0.4 to 0.435, k 2 = 0.09, 
A + = 26, ki* = 0.0168, kg = 5, kg = 0.78, and Np r g = 0.95. These values rep- 
resent the classical values generally used for equilibrium, high Reynolds number, 
turbulent flows (see ref. 12); however, it should be noted that although the 
assigned values are sufficient over a broad range of flow and wall boundary con- 
ditions, modifications to these values are required for certain classes of flow 
as discussed in previous sections of the present review. Typical turbulent 
Prandtl number variations are given in references 273 to 275 and in figure 32. 

The accuracy of the various eddy-viscosity models can only be assessed by 
careful numerical experimentation in which numerical results are compared with 
experimental data. These models have been compared over a wide range of test 
conditions for two-dimensional high Reynolds number, turbulent boundary-layer 
flows; however, their extension to general three-dimensional flows still requires 
extensive study. The models (eqs. (75), (81), and (82)) have been shown to accu- 
rately predict the data of Rainbird (ref. 276) for a sharp right circular cone 
(adiabatic wall) at angle of attack. (See refs. 271 and 277.) However, more 
demanding three-dimensional flows should be calculated and the numerical results 
compared with data. The previously mentioned problem with e x 2 must receive 
careful attention. It is not sufficient to simply assume e X 2 is some percen- 
tage of e x1 ; this has been indicated in reference 265 for infinite swept-wing 
flow. 


1 + 5.5 


x 3 ) 

^ x 3,ey 


Numerical Solution Procedures 

A number of research papers have been published over the past few years 
dealing with numerical procedures for solving the three-dimensional boundary- 
layer equations for -laminar and turbulent compressible flow. Currently, there 
is no general method for solving the three-dimensional boundary-layer equations 
(ref. 235); however, significant progress has been made in the development of 
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such procedures by Kendall et al . (ref. 245), Cebeci et al. (ref. 244), and Wort- 
man (ref. 162). A limited number of problems have been solved by applying special 
techniques to each specific problem (stagnation point flow, leading-edge attach- 
ment flow, swept infinite wings, and sharp and blunt cones). For these particu- 
lar problems it has been possible to use simple geometrical relationships (coor- 
dinate system) and similarity type transformations to reduce the complexity of 
the governing equations. As previously mentioned, two of the more difficult 
numerical problems are ( 1 ) the generation of the surface coordinate systems and 
the required metric coefficients h-| , h2, and h^ and (2) the proper treatment 

of the inviscid edge conditions. 

The properties of the three-dimensional boundary-layer equations were first 
studied in detail by Raetz (ref. 248). It was noted that the governing equa- 
tions are hyperbolic in planes parallel to the wall boundary as opposed to para- 
bolic for two-dimensional flow. Raetz introduced the concept of the influence 
principle into the three-dimensional boundary-layer literature which results 
directly from the hyperbolic character of the equations. The zones of influ- 
ence and dependence for the three-dimensional boundary-layer equations have also 
been examined in detail by Wang (ref. 256). The hyperbolic character of the 
the equations has been utilized by Bradshaw (ref. 278) through the application 
of the well-known method of characteristics, the wall region of the flow being 
patched to the characteristics solution for the outer region. 

In order to obtain the correct solution to the three-dimensional equations, 
the solution procedure must correctly account for the zone of dependence. The 
stability of the procedure and the region where the solution can be obtained 
with specified initial conditions are determined by the zone of dependence. 

This zone thus prescribes the minimum amount of initial data that must be speci- 
fied in order to advance the solution. Kitchens, Gerber, and Sedney (refs. 254 
and 257) have made detailed studies of this requirement during which they sys- 
tematically violated the zone of dependence requirement. Their results clearly 
demonstrate that errors accumulate and grow in the numerical solution if the 
zone of dependence is not satisfied. The most interesting result, presented in 
reference 254, is that satisfying the zone of dependence criteria is not suffi- 
cient to assure stability in all cases. The conclusion is drawn that the zone 
of dependence concept is not necessarily the same as the concept of stability. 

This particular result of reference 254 should be further studied since most 
investigators have treated these concepts as either the same or closely related. 

For two-dimensional boundary-layer flows the zone of dependence is automati- 
cally satisfied, the equations are parabolic, and the procedure can march down- 
stream, provided the necessary initial, boundary, and edge conditions are speci- 
fied. For three-dimensional boundary-layer flows it is necessary to internally 
adjust the mesh aspect ratio (see ref. 254), such that 

£ max A Ue A (87) 

Ax 1 \ U 1 / 

over all interior points 

Then the largest allowable Ax-| from one data plane to the next is 
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(88) 


Ax-| < a 


&X2 

max (I u 2 1 / u 1 ) over entire data plane 


where a ;£ 1 . The requirement that equations (87) and (88) be satisfied in 
the solution plane introduces numerical complexity for the three-dimensional 
boundary- layer equations. 


The numerical procedures currently being applied to solve the three- 
dimensional boundary-layer equations have resulted, in general, from a direct 
extension of procedures developed during the 1960's for two-dimensional and axi- 
symmetric boundary-layer flows. The present review is not intended to give a 
complete survey of the extension of these methods to three-dimensional boundary- 
layer flows. Reviews of three-dimensional laminar boundary layers have been pre- 
sented by Cooke and Hall (ref. 279), Stewartson (ref. 280), Crabtree, Kiichemann, 
and Sowerby (ref. 28 1 ) , Mager (ref. 282), and Hansen (ref. 283). A recent 
review of two- and three-dimensional flows has been presented by Blottner 
(ref. 235). Solutions of the three-dimensional boundary-layer equations have 
been considered by Fannelop and Humphreys (ref. 253). The reader interested in 
the initiation of numerical research in three-dimensional boundary-layer flows 
should carefully review these references , as well as The European Research Pro- 
gramme on Viscous Flows (ref. 266) and the results of Euromech 60 (ref. 284). 

The purpose of the present section is to present three solution procedures which 
are representative of current programs for three-dimensional turbulent boundary- 
layer flows; two of these procedures are currently being developed as production 
codes for arbitrary configurations. 

Crocco variables .- The primary advantage of a Crocco-type transformation is 
that the solution domain in the normal coordinate is bounded between the defi- 
nite limits of 0 g £ ^ 1 . The procedure is also attractive in that the grid 
points can be uniformly placed in the velocity plane and still satisfy the 
demanding mesh point distribution for the wall region of turbulent flows. The 
only disadvantage of the method appears to be the restraint that velocity over- 
shoot in the ui velocity component is not allowed; that is, u-)/ui >e must 
increase monotonically from the specified wall value (slip at surface can be 
allowed) to unity at the edge boundary. The transformation has been used for 
laminar flows (see refs. 255, 258, 262, and 285) and for turbulent flows over 
cones at angle of attack by Harris and Morris (ref. 1 68 ) . A discussion of 
Crocco-type variables is presented by Davis (ref. 286). 

A three-dimensional compressible turbulent code has been developed at 
Langley Research Center for Crocco-type variables (ref. 168). The program is 
currently restricted to analytic geometry; however, work is currently under way 
to extend the code to a general curvilinear coordinate system for arbitrary 
geometry . 

In the subject program, equations (60) to (63) are first nondimensional- 
ized, and a similarity-type transformation is introduced for the x^-coordinate 
and velocity as follows: 

ho = l L ho (89) 

P 
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1 


(90) 


u 3 


L 



where is the reference velocity and £ □ \y3Tj" for a sharp cone. The 

metric coefficients h2 and h 3 are arbitrary functions of the coordinates. 
The equations are next cast into Crocco form as follows: 


t = 


u-| > 
U 1 > ey 


1/2n 


(91) 


The exponent n can be selected to minimize the number of nodal points in the 
^-direction (normal to the wall boundary). The continuity and x-|-momentum equa- 
tions are combined to form the shear equations, where the shear parameter $ is 
defined by 


* = -iT 3 £H r >) (92) 

Consequently, ui/u-| je is replaced by $ as a new dependent variable, and 
H = U3/u-| je is uncoupled from the system of equations. The system therefore 
reduces to three coupled nonlinear partial differential equations in ©, G, 
and $, together with an explicit algebraic relationship for H. The system 
assumes the following form: 

^ + a-i + apw + ap + an 

ac 2 1 3 ? 3 4 K 

where u> represents 0, G, and $, and 
linear coefficients. 


+ a 5 ^ = 0 (93) 

a-| , a 2 , a 3 , an, and 05 are non- 


The boundary conditions on equation 


(93) are as follows: 


When C = 0 

G = 1 0=1 $ = 0 

When r, = 1 


0 = 0 W or 
G = 0 



f(5,n) I 





-|a,H . a 2 f 



(94) 


(95) 


where_ a-| and a2 are functions of geometry and the inviscid edge conditions 
and l = p/T. 


Equation (93) is solved in an iterative mode with a marching implicit 
finite-difference technique suggested by Dwyer (ref. 287) and modified by Krause 
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(ref. 288). The method is second-order accurate and unconditionally stable (con- 
ditional stability for reverse cross flow; see ref. 288). For turbulent flows a 
minimum of two mesh points in the C-plane must be located in the viscous sub- 
layer; consequently, careful consideration must be given to the placement of the 
mesh points in the ^-coordinate in order to minimize computer storage and pro- 
cessing time requirements. Two approaches to the minimization problem have been 
considered: the specification of a variable mesh point distribution for n = 1 

of the form ACk+lMCk = Constant for k = 2, 3, • • ., N - 1 (geometric pro- 
gression), and the selection of n = 3/4 which appears optimum for turbulent 
flows. Variable mesh-point distributions are also used in the £- and n-planes 
to minimize the computer processing time and storage requirements. A schematic 
of the difference molecule is presented in figure 55(a). Equation (93) is writ- 
ten at the point (i- 1 / 2 ,j,k) and solved for the values of the dependent varia- 
bles 0, G, $, and H at the point (i,j,k). The difference relationships 
used in the procedure for the Krause (ref. 288) "zig-zag" scheme are presented 
in reference 27 1 . 

When a converged solution cannot be obtained at the most leeward plane, 

<|> = 180° (for example, separation on leeward surface), a cubic Crank-Nicolson 
differencing scheme (see fig. 55 (b)) is used at the maximum n-station (ref. 261 ). 
If this procedure were not incorporated into the program logic, one n-station 
would be lost for each additional AC step because of the Krause differencing 
scheme. 

The marching procedure cannot be initiated without the existence of two 
orthogonal initial data planes. For a sharp right circular cone these planes 
of initial data are generated directly from the governing equations by using a 
second-order Crank-Nicolson scheme for the two planes (£ = 0, 0 g n £ n m ax and 

0 % C £ £ max , n = 0) where similarity exists. A discussion of problems associ- 
ated with obtaining initial data planes for general configurations is presented 
in reference 260. 

Substitution of the difference quotients (see ref. 271) into equation (93) 
results in a system of coupled algebraic equations whose coefficient matrix is 
of tridiagonal form which can be efficiently solved for the dependent variables 
(Thomas' algorithm). The primary problem associated with equation (93) is that 
the coefficients (a-| , a 2 , etc.) are highly nonlinear. The shear equation con- 
trols the convergence rate of the numerical procedure (the number of iterations 
required as the system is sequentially iterated). Equation (93) can be written 
for $ as 


+ Am + £l\w - liI$ + ^2 + ^39^ + 543i-0 (96) 

a c 2 \z $ $2 dK $2 9n 

where the coefficients B i , 82* 83 , and 84 are functions of geometry, invis- 

cid edge conditions, and previous iterate values of the dependent variables 
u 1 /ia 1 , e and © and their derivatives. The problem is further complicated by 
the inclusion of the turbulence models (eqs. (75) to (85)) since in the trans- 
formed plane $ appears explicitly in the transformed relationships. Conse- 
quently, the coefficients ( B 1 , 82 * etc.) also depend on $ for turbulent flows 

(for laminar flows this dependence is removed). The system of equations will 
not converge if the shear equation is written as shown in equation ( 96 ) because 
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of the term. Convergence can be achieved by using a Taylor’s series expan- 
sion of about the previous iterate value $q; that is, 

2 = ±_(2 - *_\ + 0 ($ 2 ) 

$ $G \ / 

Substitution of equation (97) into equation (96) yields 

9^$ + Am + + /_ - 5i_U 4. 2B 2 ^3 5^ - n 

a; 2 $2 )K \ ;2 $ G 2j $ G $2 3; $2 9n 

Equation (98) converges in an average of five to seven iterations for 
high Reynolds number turbulent flow. The wall boundary condition on $ (see 
eq. (95)) also presents a problem since $ w is unknown; however, the wall deriv- 
ative relationship can be directly incorporated in the iterative solution proce- 
dure. In principle, it should be possible to reduce the average number of itera- 
tions subtantially to a maximum of three. Research continues in the areas of 
(1) restructuring equation (98), (2) treatment of the $ w wall boundary condi- 
tions, and (3) the problems associated with $ in the transformed turbulence 
models. Note however that the procedure requires essentially the same process- 
ing time per mesh point (0.002 sec) as the Cebeci-Keller box method (refs. 289 
to 291) and that this time may be substantially reduced through convergence 
accelerator procedures and/or the inclusion of Newton-Raphson iteration. 

The numerical procedure and turbulence models have been applied to a number 
of flows (current geometry limited to sharp right circular and elliptic cones). 
Numerical results compared with experimental wall and profile data for a cone 
(ref. 276) with a 12.5° semiapex angle at an angle of attack of 15.75° are shown 
in figure 56. The free-stream Mach number, total pressure, and total temperature 
for the experiments were 1.8, 172.4 kN/m 2 , and 294 K, respectively. _ Transition 
was assumed to be initiated and completed in the region 0.03 ^ x-]/l £ 0.08 
(L = 105.6 cm). The adiabatic wall boundary condition was imposed on the energy 
equation (see eq. (95)); that is, (9©/9c;) w = 0. No experimental data were 
input into the viscous flow solution. The inviscid pressure distribution 
p e = f(£,ti) was obtained from a numerical solution of the three-dimensional 
inviscid flow equations. 

The numerical results for u-|/u-| e , G, and © are compared with experimen- 
tal data in figure 56 for circumferential locations of if, 0° and 135°. In 
order to evaluate the effect of nodal-point spacing in the £-plane, a parametric 
study was made for N = 301, 201, 101, 61, and 21 with ACk+l^^k = 1-02. The 
results for N = 301 and 201 were essentially identical, and those for N = 101 
were within 0.5 percent of the N = 301 results. The agreement between the 
experimental and numerical results is very good for 301 points and, in general, 
good for as few as 21 points. The two turbulence models (eqs. (75) to (85)) 
produced essentially identical results. The two-layer model results presented 
in figure 56 are for N = 301; however, the two-layer results for N = 21 were 
essentially identical to the N = 21 results of the single-layer model. 

The results for A^ k+ -|/ACk = I- 02 presented in figure 56 were obtained for 
n n 1. (See eq. (91).) Numerical results for A£ k +l/A£k = 1 (equal mesh dis- 
tribution), N = 21, and n = 3/4 were in slightly better agreement with the 


(97) 


(98) 
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N = 301 results for the variable mesh point distribution; however, the differ- 
ence between the two cases for N = 21 was minor. The major points that 
should be noted in these comparisons are that the numerical procedure is effi- 
cient and accurate and that the turbulence models are satisfactory for high 
Reynolds number equilibrium turbulent boundary-layer flows. The Crocco-type 
transformation and the numerical procedure allow the generation of accurate solu- 
tions for a minimum of 21 points normal to the wall boundary. The computer code 
requires 6OOOO3 storage (the i-1,j,k data plane is stored on disk) and approx- 
imately 0.002 second per grid point processing time on a CDC-6600 computer sys- 
tem. Current studies indicate that it may be possible to substantially reduce 
the processing time through convergence accelerators for the shear equation 
(eq. (98)) and/or the inclusion of a Newton-Raphson iteration procedure. The 
current program is comparable in both storage and processing time with the 
Cebeci-Keller box method (ref. 260). 

Arbitrary wings .- An accurate and efficient computer code for the three- 
dimensional boundary- layer flow over wings is required for the design and evalu- 
ation of supercritical wings and laminar flow control surfaces (U3 w = g(x-|,X2)). 
Studies of wing geometry specifications have indicated that a nonorthogonal sur- 
face coordinate system is optimum from the viewpoint of the design engineer. 
Cebeci et al. (ref. 260) have developed an efficient and accurate procedure for 
solving the three-dimensional boundary-layer equations for laminar, transi- 
tional, and turbulent perfect-gas flow over general wings. The advantage of 
this geometry routine are as follows: (1) Calculation of the coordinate system 

(5,ti»hi, . . . ) for each angle-of -attack or flow condition change is eliminated; 

(2) discontinuities associated with patched coordinate systems are eliminated; 

(3) the method is optimum from the viewpoint of the design engineer. (See 
sketch (i).) Only one disadvantage is encountered in that the nonorthogonal 
metric tensor results in additional terms appearing in the system of equations. 
However, this increase in terms is minimal if cross-flow diffusion terms are 
neglected. Preliminary results obtained in reference 260 for a typical super- 
critical wing indicate computation times of approximately 30 sec on an 

IBM 370/165 computer system for a 30 x 20 x 20 grid. 
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The governing equations are written for a nonorthogonal surface coordinate 
system for ease of placement of the mesh distribution over the wing surface. The 
wing is defined in an orthogonal x-|,X2,X3 coordinate system where x-j is in 
the direction of the airplane's longitudinal axis, X2 is in the spanwise direc- 
tion, and X3 is orthogonal to the plane of x-| and X2- The wing is described 
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by a series of airfoil sections in planes of X2 = Constant. The equations 
relating the airfoil section specification and the x^-coordinate system are 
presented in reference 260 in terms of percentage chord. The equations required 
for the calculation of the nonorthogonal surface mesh and the procedure with 
which they are solved are also presented in reference 260. 


The x^-momentum equation is as follows (see eq. (61) for orthogonal system): 


pu-j 3u-| 
hi 3xi 


+ puo - pKiui 2 ctn 0 + pK2Ug 2 esc 0 + PK12U1U2 
h2 3x2 ^ x 3 

_ _ esc 2 0 3p . ctn 0 esc 0 3p . 3 ( u 9u 1 _ ou i u ’\ 

hi 3xi h 2 9x 2 3 x 3\ 9x 3 7 


(99) 


where 0 represents the angle between the coordinate lines xi = Constant and 
X2 = Constant. For an orthogonal coordinate system, 0 = ir/2 and equation (99) 
reduces identically to equation (61) for h^ = 1 • The geodesic curvatures are 
given by the relationships 


Ki 


1 

hih2 sin 0 


-. 9 -(h? cos 
3xi 


0) 


3hi 

3x2 


( 100 ) 


K 2 


B 


1 

h^2 sin 0 


_ 9 — (hi cos 
3x2 


0 ) 


3h2 

3xi 


( 101 ) 


The parameter K12 is defined as follows: 


K12 


1 

sin 0 



+ 1_ \ + cos 0 ( k ? 

hi 3xij \ 


J_ 30 \ 
h 2 3x2/ 


( 102 ) 


It should be noted that the addition of correction terms for nonorthogonal coor- 
dinates to the system of equations is minor for the classical three-dimensional 
boundary- layer equations; however, if one wants to include the cross-flow diffu- 
sion terms in the X2-momentum equation, then a significant number of additional 
correction terms would be required. (See ref. 12.) 


The system of governing equations is transformed through the introduction 
of a two-component vector potential and a similarity-type transformation. (See 
ref. 260.) The Cebeci-Keller box method (refs. 162, 289, and 290) is used to 
solve the resulting system of equations. One of the basic concepts of the box 
procedure is to rewrite the system of equations as a first-order system of 
partial-differential equations. Consequently, derivatives of some quantities 
with respect to the x^-coordinate must be introduced as new unknown variables. 
Derivatives with respect to all other variables occur only to the first order 
because of the boundary-layer approximations. Centered difference quotients and 
averages at the midpoints of net rectangular and net segments are used in order 
to produce second-order-accurate finite-difference equations. The method is 
unconditionally stable for positive cross flow but appears to be unconditionally 
unstable for negative cross flow. The equations are highly nonlinear and impli- 
cit in structure. Newton's method is used to solve the system; a block- 
tridiagonal factorization scheme is used which is efficient and stable. The 
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numerical formulation of the system of equations is presented in detail in ref- 
erence 28. Numerical results are compared with the data of East and Hoxy 
(ref. 291) in figure 57 for several x-j/L stations. Numerical results for a 
supercritical wing calculation are presented in figure 58. 

Reacting gas flows .- The Aerotherm Division of the Acurex Corporation is 
currently developing a laminar, transitional, and turbulent three-dimensional 
boundary-layer computer code for application to arbitrary configurations such as 
the space shuttle. (See ref. 245.) The numerical procedure is an extension of 
the method presented in reference 238. This particular code includes equilib- 
rium and/or frozen flow chemistry. A body-oriented orthogonal surface coordi- 
nate system is used to describe the vehicle surface. The computer code is com- 
patible with the accurate three-dimensional supersonic inviscid flow field 
program developed by Marconi, Yaeger, and Hamilton (ref. 246). Three- 
dimensional entropy- layer effects are included in the solutions. This is of 
particular importance for vehicles such as the space shuttle. 


The governing equations are numerically solved in primitive variables with 
suitable stretching and normalization in the x^-coordinate . Similarity trans- 
formations are not used since the program will in general be applied to highly 
nonsimilar flows. The governing equations, in nondimensional form, are pre- 
sented in reference 245. The solution domain is covered by a mesh point network 
shown in figure 59. The initial data planes are assumed known at x = x^; the 
solution is to be obtained at X£ + -j = = Ax^. The functional form of the 

derivatives is specified and substituted into the governing equations together 
with the required boundary conditions. This procedure results in a system of 
algebraic relations between the unknown dependent variables at the nodal point 
in the plane located at x^ + -j . Normal to the wall boundary, the dependent vari- 
ables are represented by a splined Taylor series between each mesh point as fol- 
lows (where the prime denotes differentiation with respect to X3): 


„ . A « AXO 2 

u j+1 = u j - u j Ax 3 + u j 


u J+1 = uj + Uj Ax 3 


(103) 


Where j represents the grid index (j = 1, 2, . . ., JMAX) . The equations are 
integrated with respect to X3 between each mesh point. This insures that the 
conservation laws are satisfied exactly between the mesh points and greatly sim- 
plifies the calculation of the diffusion terms since all second derivatives are 
eliminated. 


In the cross-flow direction (x2~coordinate) second-order-accurate centered 
difference quotients are used to replace the cross-flow derivatives; that is, 
for equal nodal spacing 

= u i, J,k+1 - u i, j,k-1 ( 104 ) 

W’A.j.k Ax 2,K*’ 

where k is the index for the X2 mesh point distribution. Axial derivatives 
are approximated by the backward difference quotients of the form 
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(105) 



u i ± 1 i j,k " u i,j,k 


Ax 


1,i 


where i is the index associated with the mesh point distribution in the 
xi -coordinate. 


A fully implicit solution is obtained in the procedure. Equations ( 103 ) 
and (104) are evaluated in the x-| ^ + -|-plane; introduction of the boundary con- 
ditions results in a nonlinear system of algebraic equations for the primary 
variables (u, p, etc.) and secondary variables (u', p' , etc.) at each nodal 
point in the x-| ^ + -|-plane. The resulting system of algebraic equations is 
solved by the Newton-Raphson procedure. The solution reduces to the problem of 
inverting a large matrix and thereby solving the linear system of algebraic equa- 
tions for the unknowns. The large matrix is treated as a system of submatrices. 
The matrix is block tridiagonal in structure since cross-flow derivatives are 
represented by centered difference quotients. The submatrices along the main 
diagonal are dense, whereas the off-diagonal submatrices are sparse. A schema- 
tic of a typical diagonal submatrix for one x-| -plane is presented in figure 4 
of reference 245. 


Data required for the initial planes are obtained directly in the solution 
procedures at stagnation points and leading-edge attachment lines. Similarity 
variables are used for the generation of these solutions. The entropy layer 
effects and the procedure through which they are included in the numerical pro- 
cedure are discussed in detail by Kendall et al. (ref. 245). The eddy-viscosity 
model currently used is the single-layer model (see eq. (75)), where the Lewis 
and Prandtl numbers are assumed constant. More general models will be incorpo- 
rated into the code once it is completely operational and detailed studies have 
been made of comparisons of numerical results with experimental data. 

A comparison of numerical results with experimental data (ref. 292) for 
flow over a sharp cone at angle of attack is presented in figure 60. Solutions 
for the flow over a space shuttle are currently being obtained. 


Status of Three-Dimensional Boundary-Layer Computational Techniques 

Significant progress has been made in the development of numerical tech- 
niques for solving the compressible three-dimensional turbulent boundary- layer 
equations (see appendix, table A2); however, more efficient and reliable methods 
must be developed for application to general aerospace vehicles. This becomes 
clearly apparent when one considers the prohibitive computer processing times 
which can occur for either complex reacting gas flows or turbulent flows over 
large-scale aerospace vehicles (ref. 260). Early numerical research in the area 
of compressible, three-dimensional boundary- layer flows was primarily directed 
towards developing second-order-accurate, stable numerical schemes for solving 
the laminar equations for flows where simplifying assumptions were possible. 
Current research is directed towards ( 1 ) improving these numerical procedures 
for application to complex flows, (2) developing more general coordinate sys- 
tems, and (3) detailed studies of turbulence modeling procedures for three- 
dimensional turbulent-boundary- layer flows. Reynolds stress models require accu- 
rate experimental data where the stress components are measured. The wealth of 
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two-dimensional experimental data which made possible the successful development 
of two-dimensional eddy-viscosity models may provide little reliable guidance 
for complex three-dimensional boundary- layer flows. 

The problems associated with solving the three-dimensional inviscid equa- 
tions for complex vehicles are rapidly being solved; however, efficient and reli- 
able procedures for automatically coupling the inviscid and boundary- layer codes 
must be developed in order to avoid time-consuming and difficult data manipula- 
tion. Compatible coordinate systems must be developed for the inviscid and vis- 
cous codes in order to reduce or completely eliminate the inherent errors associ- 
ated with data interpolation and smoothing between coordinate systems. Current 
studies and experience indicate that a nonorthogonal curvilinear coordinate sys- 
tem is optimum for general three-dimensional boundary-layer flows. 

Flexibility must also be programed into the boundary-layer codes which will 
assure that the optimum difference scheme is utilized at various locations in 
the solution domain in order to automatically satisfy the zone of dependence 
requirements as well as to maximize the region over which the solution may be 
obtained. Care must be exercised in order to determine whether the classical 
three-dimensional boundary- layer equations are valid in a specific area where 
the solution is required; that is, boundary region flow, inflow regions, and 
regions of separated flow require special treatment. For these cases either the 
parabolized Navier-Stokes or the full Navier-Stokes equations must be solved. 

A number of difficult problems remain to be solved before compressible, 
three-dimensional, turbulent, boundary- layer codes are developed to the current 
confidence level of existing two-dimensional codes. An optimistic review of cur- 
rent research programs, together with the increasing capabilities of digital com- 
puter systems and the maturing of three-dimensional inviscid flow field codes, 
indicates that general purpose, three-dimensional boundary-layer codes for com- 
pressible turbulent flows will become available in the near future. 


RESUME 

Compressibility Influence on Turbulent Boundary-Layer Shear Stress 

From comparisons of high-speed data with low-speed closure procedures using 
variable mean density, there does not appear to be any appreciable influence of 
compressibility upon turbulent shear stress modeling in compressible turbulent 
boundary layers, even for extreme cases such as Mach 14 to 20 with a change in 
density across the layer of up to a factor of 100. Other evidence of an appar- 
ent lack of any compressibility-caused new physics which may alter the shear 
stress for the compressible boundary-layer case includes the following: 

(1) Fluctuation Mach number is generally less than 1; (2) the shear stress dis- 
tribution through the boundary layer is not a function of Mach number for zero 
pressure gradient flows; (3) the Morkovin hypothesis is valid up to Mach 5 
(based on fluctuation data); (4) profile N power is not a function of Mach num- 
ber, at least up to Mach 10; (5) the nondimensional burst period is approxi- 
mately the same as that for low speed. 
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Computational Capability of Existing Procedures 

Compressible equilibrium and near equilibrium boundary-layer flows, at 
least up to Mach 20, can be computed fairly accurately using mean field methods, 
provided that the influence of such items as low Reynolds number effects, pres- 
sure gradient, and wall blowing are properly accounted for through adjustments 
in the modeling constants. By and large, these adjustments can be obtained from 
low-speed data. Several calculation methods are available for the nonequilib- 
rium case, but these have not yet received sufficient exercise on compressible 
highly nonequilibrium flows (primarily as a result of lack of data) to determine 
their relative or absolute accuracy. 


Important Unresolved Issues 

The following is a brief listing of what are the more important unresolved 
issues (or research frontiers) in the calculation of compressible turbulent 
boundary layers, as discussed in the present paper: 

1 . Considerable further development and calibration of nonequilibrium meth- 
ods is required for compressible flows. Further experimental data in highly non- 
equilibrium flows must be obtained before this effort can be really meaningful. 

2. Definitive experiments are needed to determine directly the influence of 
compressibility upon shear stress production in compressible turbulent boundary 
layers (conditional sampling measurements at high Mach numbers and measurements 
of fluctuating pressure terms in Reynolds stress equations). 

3. Further definition is needed of the boundaries between boundary-layer 
and nonboundary-layer behavior; that is, where can boundary-layer methods be 
reasonably expected to work? Bradshaw addresses this with his "extra rates of 
strain," but definitive guidelines are needed (such as how much concave longitu- 
dinal curvature is necessary before longitudinal vortices develop and the flow 
is no longer a boundary layer). Obviously one can employ Navier-Stokes codes 
and complete Reynolds stress modeling to attack such situations, but currently 
the boundary-layer codes are so much faster to run on computers (as compared 
with the Navier-Stokes solvers) that they should be used whenever possible; this 
necessitates having well-known limits within which accurate answers can be 
obtained (without new physics making either the basic equation set or the turbu- 
lent modeling inapplicable). 

4. Obviously, three-dimensional compressible turbulent-boundary-layer calcu- 
lations are in an early stage, and considerable experimentation is necessary 
before the efficacy of mean field methods for the three-dimensional case is 
determined. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
April 7, 1977 
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APPENDIX 


SUMMARY OF CALCULATION PROCEDURES FOR NONSIMILAR TWO- AND 
THREE-DIMENSIONAL COMPRESSIBLE TURBULENT BOUNDARY 
LAYERS (FINITE DIFFERENCE, FINITE ELEMENT, AND 
WEIGHTED RESIDUAL METHODS) 

The purpose of this summary (or catalog) of numerical prediction methods is 
to indicate the wide variety of such procedures which are available. (Table A1 
has 31 entries.) In most cases information for the various entries were 
obtained from the individual authors, who kindly filled out and returned data 
sheets on their methods. Several authors have made a special effort to optimize 
the numerical solution procedure (e.g., refs. 9, 12, 162, and 1 63 ) to reduce 
the required machine time and storage per case. Most of the procedures have 
detailed user’s manuals, and in many cases the codes are available either with- 
out charge or for the cost of mailing. Many of the special effects treated by 
the various methods (such as nonequilibrium or equilibrium chemistry, transi- 
tion, roughness, etc.) are indicated in tables A1 and A2, as is the fundamental 
closure approach. 

If these tables serve no other function, they should at least cause 
researchers to think quite carefully before producing the 32nd entry; that is, 
several of the procedures are quite similar, and future research should obvi- 
ously concentrate on using the best numerical methods to investigate such items 
as nonequilibrium flows, rather than on producing another compressible turbulent- 
boundary-layer deck. 
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TURBULENT BOUNDARY LAYERS 


APPENDIX 
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TABLE A1 . - Continued 
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Figure 1 Typical distribution of p with y at high Mach number 




Figure 3.- Compa rison of T 3tatio with u ,2 /2c p at high Mach numbers. 
VU'^/Ug 2 = 0.05; T w /T t «* 1; y = 1.4. 
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Figure H.- Illustration of disturbances upstream of transition. 
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M = 13.6 
e 



GROWTH OF TURBULENCE IN A HYPERSONIC BOUNDARY LAYER BASED ON 
SCHLIEREN PHOTOGRAPHS AND SURFACE HEAT TRANSFER 

Figure 5.- Flow field schematic of precursor transition for Mach 18 cone 

experiments (taken from ref. 57). 
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Figure 6.- Typical velocity profile at beginning of conventional transition 
region for M » 1 (from ref. 58). M e = 8.96; T w /T t = 0.656; 
x tr = 54.4 cm; end of transition at approximately 109.3 cm. 
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Figure 7.- Influence of precursor transition effect 
on wall heating (from ref. 60). 



Figure 8.- Height of critical layer as function 
of Mach number (from ref. 62). 
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Figure 9.- Influence of Mach number on fluctuating wall pressure levels 
for turbulent boundary-layer flow (from ref. -47 ) . 



Figure 10.- Distribution of fluctuating static pressure through a turbulent 
boundary layer in high Mach number flow. M e = 9.4; R e,0 = 36 800; 

T w /T t h 0.38. (Reprinted from ref. 62 with permission from Cambridge 
University Press.) 
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M ~ 20 



Figure 11.- Illustration of difference between pitot (or density) and 
velocity boundary- layer thicknesses at high Mach number (from 
ref. 47). density is y-location where p/p e = 0.995 and ^velocity 
is y-location where u/u e = 0.995. 
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Figure 12.- Prediction of flat-plate skin friction by transformation methods 
of Baronti and Libby (ref. 75) and Economos and Boccio (ref. 78). 
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Figure 13.- Prediction of flat-plate skin friction and heating 
by transformation method of Coles (ref. 15). 
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(a) Economos and Bocoio (ref. 78). 
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(b) Lewis et al. (ref. 77). 


Figure 14.- Prediction of boundary-layer characteristics on a body 
of revolution (data from ref. 83 ). 
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Figure 15.- Prediction of flat-plate skin friction and heating by. correlation 
method of Van Driest (ref. 87). M = 0 to 10; T w /T^ = O.lM to 0.7. 
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Figure 16.- Prediction of flat-plate skin friction and heating by correlation 
method of Spalding and Chi (ref. 81). M = 0 to 10; T w /T^ = 0.1 4 to 0.7. 
(R.A.F. , Reynolds analogy factor.) 
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HEAT TRANSFER 

(VON KARMAN'S R.A.F. ASSUMED) 



Figure 17.- Prediction of flat-plate skin friction and heating by correlation 
method of Eckert (ref. 72). M = 0 to 10; T w /T t = 0.14 to 0.7. (R.A.F., 

Reynolds analogy factor.) 



Figure 18.- Prediction of turbulent skin friction on conventional flight 
vehicles by Eckert's method (ref. 72). 
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Figure 19 •- Prediction of turbulent heating on sharp and blunt tipped 
axisymmetric bodies in flight by Eckert's method (ref. 72); data 
from reference 96 . 
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(b) Waisted body of revolution (experimental data from ref. 83). 


Figure 20.- Prediction of boundary layer properties for turbulent flows with 
pressure gradient by Flaherty's integral method (ref. 124). 
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(a) Waisted body of revolution (experimental data from ref. 83). 



(b) Flat cylinder (experimental data from ref. 139). 


Figure 21 .- Prediction of boundary-layer properties for flows with pressure 
gradients by Green's integral method (ref. 135). 
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Figure 22.- Prediction of boundary-layer development on a concave surface 
by Pinckney's integral method (ref. 133). Experimental data from 
reference 141 . 
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(a) Flat plate with favorable pressure gradient (experimental 

data from ref. 138). 
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(b) Waisted body of revolution (experimental data from ref. 83). 

Figure 23.- Prediction of boundary-layer properties for flows with pressure 
gradient using Reeves' integral method (ref. 30, pp. 6-1 - 6-A2-2). 
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Ay 

Figure 24.- Effect of grid spacing on computational error (from ref. 169). 
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Figure 25.- Effect of variable grid spacing on calculated momentum thickness 
and skin friction. (Reprinted from ref. 170 with permission from the 
American Institute of Aeronautics and Astronautics. 
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10 1 10 2 6 + io 3 10 4 


(a) Flow over plates, cones, and cylinders. 



6 + 


(b) Nozzle wall flow. 

Figure 29.- Variation of U/S)maY with 6 + for various types of flow 

(frSref. 189). 
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Figure 30.- Illustration of increasing importance of low Reynolds 
number effects at high Mach number. 



Figure 31.- Effect of wall injection and positive pressure gradient 

on ( i / 6 ) max (from ref. 53). 
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Figure 32.- Variation of turbulent Prandtl number with y + . 



0 .2 .4 .6 .8 1.0 

y/6 


Figure 33.- Comparison of static with total turbulent Prandtl number at 
Mach 7.2. (Reprinted from ref. 177 with permission from the American 
Institute of Aeronautics and Astronautics.) 
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6 + 

Figure 3M.- Variation of Reynolds analogy factor with 6 + . 



Figure 35.- Theoretical variation of turbulent Prandtl number 
with e/u. (Reprinted with permission from ref. 12.) 
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Figure 39.- Development of total temperature-velocity profile 
along a nozzle wall (data from ref. 212). 



max 


= ((S' 1 ) 


CALCULATIONS 


Figure 40.- Influence of wall blowing on turbulent skin friction 

(from ref. 217) . 
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LOCAL MACH NUMBER 

Figure 41.- Extent of transition zone in high speed flow. (Reprinted from 
ref. 219 with permission from the American Institute of Aeronautics and 
Astronautics. ) 
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Figure 42.- Effect of Mach number on streamwise intermittency . (Reprinted 
from ref. 59 with permission from the American Institute of Aeronautics 
and Astronautics.) 
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Figure 43.- Comparison of transitional calculations including precursor and 
low Reynolds number effects with data (from ref. 58). The profile is 
near center of transition region (x = 78.7 cm). 
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Figure 44.- Effects of variable entropy on wall heating data and calculations 
(from ref. 221). M = 10; T t = 1111 K; T w = 294 K; A+ =26; K = 0.435; 
( l /6 )may = 0 - 09 ; Npp^ = 0.9* 
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Figure 46.- Flow conditions for data in figure 45. 
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Figure 49.- Effect of fluctuating stream velocity on a-j (from ref. 233) 
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Figure 50.- Values of a-| through a shock — boundary-layer interaction 

at Mach 3.88 (from ref. 234). 




Figure 51 Calculation of u|uj with a turbulent boundary layer using 
nonequilibrium turbulence modeling (from ref. 156). 



Figure 52.- Calculation of u'v' with a turbulent boundary layer using 
nonequilibrium turbulence modeling (from ref. 156). 
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Figure 56.- Comparison of three-dimensional prediction scheme with data (from 
ref. 276) for a cone at angle of attack in a supersonic stream (d = 2.54 cm). 
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BOUNDARY LAYER IS COVERED BY A NODAL NETWORK 



Figure 59.- Mesh point network for integral-matrix procedure. 




Figure 60.- Comparison of numerical results with surface heating data (from 
ref. 292) for a cone at angle of attack in hypersonic flow. Mo, = 7.95; 

Y = 7.95; a = 4°. 
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